Wall CB, Ricci CA, Wen AD, Ledbetter BE, Klinger DE, Mydlarz LD, Gates RD, Putnam HM. Repeat bleaching drives coral physiotypes by environmental legacies and cellular memory. (submitted Nature Climate Change)

Background

Figure 1.. Multivariate analyses of organisms and their response to environmental change and a conceptual diagram of coral physiology and immunity responses during heat stress.

Thermal stress is challenging corals and their capacity to cope with and survive recurrent belaching events. These events may shift corals to new physiotypes (i.e., multivariate ‘trait space’) as a result of shifts in host or symbiont biology and constitutive immunity to maintain homeostasis. Here we collected data from Kāne’ohe Bay, O’ahu, Hawai’i extending across two archipelagic bleaching events and the subsequent recovery periods (four Periods). Data periods represent:
- first bleaching (October 2014)
- first recovery (February 2015)
- second bleaching (October 2015)
- second recovery (February 2016)

Figure 2 (partial). A bleached Montipora capitata coral. (PC: CB Wall)

Physiological and immunological parameters were collected from from Montipora capitata collected from (two Sites):
- Lilipuna – a shallow fringing reef west of the Hawai’i Insititute of Marine Biology
- Reef 14 – a shallow inshore patch reef in central Kāne’ohe Bay

Lilipuna has been categorized as experiencing near shore impacts from freshwater inundation, subteranean groundwater discharge, and seawater biogeochemistry influeces by long residence times, low flow, and riverine/terrigenous nutrient inputs. Additionally, pCO2 flux at this site is “moderate”, with on average ~200ppm pCO2 +/- flux

Reef 14 has seawater with a shorter residence time, little terrestrial impact, but has a greater pCO2 flux due to reef metabolism. This flux can be substantial, being > 600 ppm pCO2 in a 24h period and reaching maximum values of >900ppm pCO2.

Environmental data

Temperature

Temperature data was collected at two scales. Each of the sites had loggers placed on the reef where data was continuously recorded over 3 years. Gaps in this data exist due to logger failure or loss, however, temperatures correspond well to NOAA station data. The long term data from NOAA gives indications of the rise and fall of temperautures throughout Kāne’ohe Bay during the 2024-2016 period.
Site loggers present at each of the two reef sites.

#setwd("data/environmental/Temp data")

# Import temperature data

#2014 Oct - 2015 Jun
Lil.temp_log1 <- rbind(
        read.csv("data/environmental/Temp data/Lilipuna/SN_10487932/lilipuna_oct-dec 2014.csv"),
        read.csv("data/environmental/Temp data/Lilipuna/SN_10487932/lilipuna_dec-feb 2015.csv"))

Lil.temp_log2 <- rbind(
        read.csv("data/environmental/Temp data/Lilipuna/SN_10487932/lilipuna_feb-apr 2015.csv"),
        read.csv("data/environmental/Temp data/Lilipuna/SN_10487932/lilipuna_apr-jun 2015.csv"))

#2015 Sep - 2016 Jun
Lil.temp_log4 <- rbind(
        read.csv("data/environmental/Temp data/Lilipuna/SN_10084646/lilipuna_sep-nov 2015.csv"),
        read.csv("data/environmental/Temp data/Lilipuna/SN_10084646/lilipuna_nov-jan 2016.csv"),
        read.csv("data/environmental/Temp data/Lilipuna/SN_10084646/lilipuna_jan-mar 2016.csv"))

Lil.temp_log5<- rbind(
        read.csv("data/environmental/Temp data/Lilipuna/SN_10084646/lilipuna_apr-jun 2016.csv"))
                   

#2014 Oct - 2015 Jun
Rf14.temp_log1 <- rbind(
  read.csv("data/environmental/Temp data/Reef 14/SN_10487931/reef14_oct-dec 2014.csv"),
  read.csv("data/environmental/Temp data/Reef 14/SN_10487931/reef14_dec-feb 2015.csv"))

Rf14.temp_log2 <-rbind( 
  read.csv("data/environmental/Temp data/Reef 14/SN_10487931/reef14_feb-apr 2015.csv"),
  read.csv("data/environmental/Temp data/Reef 14/SN_10487931/reef14_apr-jun 2015.csv"))

Rf14.temp_log3 <- 
  read.csv("data/environmental/Temp data/Reef 14/SN_9893760/reef14_jun-sep 2015.csv")

Rf14.temp_log4 <- rbind(
  read.csv("data/environmental/Temp data/Reef 14/SN_9893760/reef14_sep-nov 2015.csv"),
  read.csv("data/environmental/Temp data/Reef 14/SN_9893760/reef14_sep-nov 2015.csv"),
  read.csv("data/environmental/Temp data/Reef 14/SN_9893760/reef14_nov-jan 2016.csv"),
  read.csv("data/environmental/Temp data/Reef 14/SN_9893760/reef14_jan-mar 2016.csv"))
  
Rf14.temp_log5<- rbind(
  read.csv("data/environmental/Temp data/Reef 14/SN_9893760/reef14_apr-jun 2016.csv"))
  
# set date to date class
Lil.temp_log1$Date <- as.Date(Lil.temp_log1$Date, format="%Y/%m/%e")
Lil.temp_log2$Date <- as.Date(Lil.temp_log2$Date, format="%Y/%m/%e")
Lil.temp_log4$Date <- as.Date(Lil.temp_log4$Date, format="%Y/%m/%e")
Lil.temp_log5$Date <- as.Date(Lil.temp_log5$Date, format="%Y/%m/%e")

Rf14.temp_log1$Date <- as.Date(Rf14.temp_log1$Date, format="%Y/%m/%e")
Rf14.temp_log2$Date <- as.Date(Rf14.temp_log2$Date, format="%Y/%m/%e")
Rf14.temp_log3$Date <- as.Date(Rf14.temp_log3$Date, format="%Y/%m/%e")
Rf14.temp_log4$Date <- as.Date(Rf14.temp_log4$Date, format="%Y/%m/%e")
Rf14.temp_log5$Date <- as.Date(Rf14.temp_log5$Date, format="%Y/%m/%e")

##########################

##Calibrate temperature

# Lil.temp_log1: SN_10487932   y = 1.0066x - 0.0085
# Lil.temp_log2: SN_10084646   y = 0.9978x + 0.0815
# Rf14.temp_log1: SN_10487931   y = 1.0047x + 0.0722
# Rf14.temp_log2: SN_9893760   y = 1.0044x - 0.0648

Lil.temp_log1$Temp_Calib<-(1.0066*(Lil.temp_log1$Temp)-0.0085)
Lil.temp_log2$Temp_Calib<-(1.0066*(Lil.temp_log2$Temp)-0.0085)
Lil.temp_log4$Temp_Calib<-(0.9978*(Lil.temp_log4$Temp)-0.0815)
Lil.temp_log5$Temp_Calib<-(0.9978*(Lil.temp_log5$Temp)-0.0815)

Rf14.temp_log1$Temp_Calib<-(1.0047*(Rf14.temp_log1$Temp)-0.0722)
Rf14.temp_log2$Temp_Calib<-(1.0047*(Rf14.temp_log2$Temp)-0.0722)
Rf14.temp_log3$Temp_Calib<-(1.0044*(Rf14.temp_log3$Temp)-0.0648)
Rf14.temp_log4$Temp_Calib<-(1.0044*(Rf14.temp_log4$Temp)-0.0648)
Rf14.temp_log5$Temp_Calib<-(1.0044*(Rf14.temp_log5$Temp)-0.0648)

#compile calibrated data from all time points into single files 

TIMESPAN<-rbind(Rf14.temp_log1, Rf14.temp_log2, Rf14.temp_log3, Rf14.temp_log4, Rf14.temp_log5) #all dates from period

Lil.Temp1 <-Lil.temp_log1 #  Lilipuna temp data start 
Lil.Temp2 <-Lil.temp_log2 # 2 week break in data, Feb until loggers lost
#Lil.Temp3 ==loggers lost
Lil.Temp4 <-Lil.temp_log4 # all data until March 2016 gap
Lil.Temp5 <-Lil.temp_log5 # 2016 March onward

Rf14.Temp1 <- Rf14.temp_log1 
Rf14.Temp2 <- Rf14.temp_log2 
Rf14.Temp3 <- Rf14.temp_log3 
Rf14.Temp4 <- Rf14.temp_log4                       
Rf14.Temp5 <- Rf14.temp_log5                       

#############################
## complete data for both sites across all periods of sampling (2014-2016)
Lil.Temp<-rbind(Lil.temp_log1,Lil.temp_log2, Lil.temp_log4, Lil.temp_log5)
Rf14.Temp<-rbind(Rf14.Temp1,Rf14.Temp2, Rf14.Temp3, Rf14.Temp4, Rf14.Temp5)

#############################

# Aggregate temperature data by daily mean

TIMESPAN_split<-split(TIMESPAN, f=TIMESPAN$Date 
                     < as.Date("2014-10-10", format="%F"))

Lil.split1 <- split(Lil.Temp1, f=Lil.Temp1$Date 
                      < as.Date("2014-10-10", format="%F"))
Lil.split2 <- split(Lil.Temp2, f=Lil.Temp2$Date 
                    < as.Date("2014-10-10", format="%F"))
Lil.split4 <- split(Lil.Temp4, f=Lil.Temp4$Date 
                   < as.Date("2014-10-10", format="%F"))
Lil.split5 <-split(Lil.Temp5, f=Lil.Temp5$Date 
                     < as.Date("2014-10-10", format="%F"))

Rf14.split1 <- split(Rf14.Temp1, f=Rf14.Temp1$Date 
                       < as.Date("2014-10-10", format="%F"))
Rf14.split2 <- split(Rf14.Temp2, f=Rf14.Temp2$Date 
                     < as.Date("2014-10-10", format="%F"))
Rf14.split3 <- split(Rf14.Temp3, f=Rf14.Temp3$Date 
                    < as.Date("2014-10-10", format="%F"))
Rf14.split4 <- split(Rf14.Temp4, f=Rf14.Temp4$Date 
                     < as.Date("2014-10-10", format="%F"))
Rf14.split5 <- split(Rf14.Temp5, f=Rf14.Temp5$Date 
                     < as.Date("2014-10-10", format="%F"))

##########################
##########################
TIMESPAN_mean <- aggregate(data.frame(mean=TIMESPAN_split[[1]]$Temp_Calib), by=list(Date=TIMESPAN_split[[1]]$Date), FUN=mean)

# daily means
Lil_mean1 <- aggregate(data.frame(mean=Lil.split1[[1]]$Temp_Calib), by=list(Date=Lil.split1[[1]]$Date), FUN=mean)

Lil_mean2 <- aggregate(data.frame(mean=Lil.split2[[1]]$Temp_Calib), by=list(Date=Lil.split2[[1]]$Date), FUN=mean)

Lil_mean4 <- aggregate(data.frame(mean=Lil.split4[[1]]$Temp_Calib), by=list(Date=Lil.split4[[1]]$Date), FUN=mean)

Lil_mean5 <- aggregate(data.frame(mean=Lil.split5[[1]]$Temp_Calib), by=list(Date=Lil.split5[[1]]$Date), FUN=mean)

#######

Rf14_mean1 <- aggregate(data.frame(mean=Rf14.split1[[1]]$Temp_Calib), by=list(Date=Rf14.split1[[1]]$Date), FUN=mean)

Rf14_mean2 <- aggregate(data.frame(mean=Rf14.split2[[1]]$Temp_Calib), by=list(Date=Rf14.split2[[1]]$Date), FUN=mean)

Rf14_mean3 <- aggregate(data.frame(mean=Rf14.split3[[1]]$Temp_Calib), by=list(Date=Rf14.split3[[1]]$Date), FUN=mean)

Rf14_mean4 <- aggregate(data.frame(mean=Rf14.split4[[1]]$Temp_Calib), by=list(Date=Rf14.split4[[1]]$Date), FUN=mean)

Rf14_mean5 <- aggregate(data.frame(mean=Rf14.split5[[1]]$Temp_Calib), by=list(Date=Rf14.split5[[1]]$Date), FUN=mean)

NOAA-HIMB annual temps accessed here (https://www.ndbc.noaa.gov/station_page.php?station=mokh1)

# attach the 3 years of temp data from HIMB 
# note time in file is UTC
#setwd("data/environmental/Temp data/HIMB station")

annual.temps<-rbind(read.csv("data/environmental/Temp data/HIMB station/temp_2014.csv"), 
                    read.csv("data/environmental/Temp data/HIMB station/temp_2015.csv"), 
                    read.csv("data/environmental/Temp data/HIMB station/temp_2016.csv"))

#keep these columns
annual.temps<- annual.temps %>%
  select(c(X.YY, MM, DD, hh, mm, WTMP))

#remove errant values
annual.temps<-annual.temps[!(annual.temps$WTMP > 40),]

annual.temps$Date<-as.Date(paste(annual.temps$X.YY, annual.temps$MM, annual.temps$DD, sep="-"))
annual.temps$Time<-paste(annual.temps$hh, annual.temps$mm, sep=":")
colnames(annual.temps) <- c("Year", "Month", "Day", "hh", "mm", "Water.temp", "Date", "Time")

#keep these columns
annual.temps<- annual.temps %>%
  select(c(Date, Time, Water.temp))

# make a date-time column to correct for format in UTC
annual.temps$Date<-parse_date_time(annual.temps$Date, "Y-%m-%d") # corrects date format
annual.temps$date.time<-as.POSIXct(paste(annual.temps$Date, annual.temps$Time), 
                                   format="%Y-%m-%d %H:%M", tz="UTC")

# correct time zone to HST
annual.temps$date.time<-with_tz(annual.temps$date.time, tzone = "HST")

annual.temps<-annual.temps[!c(annual.temps$date.time<"2014-01-01 00:00:00"),] #start Jan 2014
annual.temps<-annual.temps[!c(annual.temps$date.time>="2016-03-01 00:00:00"),] #end February 2016

# correct Date and Time to reflect date-time HST
annual.temps$Date<-as.Date(annual.temps$date.time)
annual.temps$Time<-strftime(annual.temps$date.time, format="%H:%M:%S")

#reorder
annual.temps<- annual.temps %>%
  select(c(date.time, Date, Time, Water.temp))

# export to csv to use in DHW calculation
write.csv(annual.temps, "output/annual.temps.csv") 

# calculate rolling mean of temp data as 'daily mean'
annual.day.temps <- aggregate(data.frame(mean=annual.temps$Water.temp), by=list(Date=annual.temps$Date), FUN=mean)

Degree Heating Weeks (DHW)

With a 27.7 °C Summer Monthly Maximum Temp (Jokiel and Brown 2004) for Kāne’ohe Bay, total DHW for 2014 and 2015 are calculated from NOAA Molu o Lo’e buoy data.
These show a maximum sustained DHW in 2014 of 7.1 (peaked from 04 Oct - 19 Nov). DHW began accumulating 2014-08-30. In 2015 10.2 DHW were sustained (peaked 12 - 19 Nov), with DHW accumulating 2014-07-01.

########## CALCULATE DHW FOR KBAY 

# Load hourly temperature dataset (2014-2016)
NOAA.temp<-read.csv(file="output/annual.temps.csv")
df<-NOAA.temp[-1] #remove column
df$date.time<-as.POSIXct(df$date.time) # correct formatting
df<-df[!(df$date.time>= "2016-03-01 00:00:00"),] # end date

# make a column for each date.hour and generate means
df$date.hour <- strftime(df$date.time, format="%Y-%m-%d %H:00:00", breaks="hour")
hour.temp.means <- aggregate(Water.temp ~ date.hour, df, mean)

# make date sequence for hourly data to identify data gaps in NOAA data
date.seq<-as.data.frame(
  seq.POSIXt(from = as.POSIXct('2014-01-01 00:00:00'), to = as.POSIXct('2016-03-01 00:00:00'), by = 'h'))
colnames(date.seq)<-"date.hour"

# merge sequence of dates generated above with actual data from NOAA
merged.hour.temp<-as.data.frame(join_all(list(date.seq, hour.temp.means), by = "date.hour", type='full'))
write.csv(merged.hour.temp, "data/environmental/Temp data/merged.hour.temp.csv")

# make DHW a matrix, then make sure even to round for DHW
date.DHW<-merged.hour.temp$date.hour
date.DHW<-date.DHW[1:(floor(length(date.DHW)/84)*84)]

# dates to export
half.wk.dates<- date.DHW[seq(1, length(date.DHW), 84)] 

# calculate temp mean at every 84 values = 3.5d mean temp
half.wks.temp<-merged.hour.temp$Water.temp
half.wks.temp<-half.wks.temp[1:(floor(length(half.wks.temp)/84)*84)]
dim(half.wks.temp) = c(84,length(half.wks.temp)/84) #dimensions match above?
half.wks.mean<-as.data.frame(.colMeans(merged.hour.temp$Water.temp, 84, 
                                       length(merged.hour.temp$Water.temp)/84, na.rm=TRUE))

# BIND TOGETHER
half.week.temps<-cbind(half.wk.dates,half.wks.mean)
colnames(half.week.temps)<-c("half.week.date.time", "Water.temp.half.week")


# Set the "Base" temperature (i.e. the mean monthly maximum, from NOAA) # 
# 27.7? (https://www.ndbc.noaa.gov/view_climplot.php?station=mokh1&meas=st) 2008-2012 maximum summer mean and Jokiel and Brown
# from NOAA Main Hawaiian Islands (https://coralreefwatch.noaa.gov/vs/data/hawaii.txt):
# Averaged Monthly Mean (Jan-Dec): 24.3713 24.0663 24.0859 24.1832 24.7437 25.4589 26.0265 26.5946 26.9824 26.8742 26.1951 25.0945
base<-27.7


###################################################################################
## Now calculate the hotspot for all sites
# Create a dataframe including time vector xi3
KBay.hotspot <- data.frame(half.wk.dates)

# Create "hotspot" for ease of doing multiple sites
hotspot<-KBay.hotspot$hotspot

# Assign hotspot value, this is (halfweekly temperature) - (base temperature)
hotspot<-half.week.temps$Water.temp.half.week-base

# Only keep hotspot temperature if the hotspot is >1 (for calculation of DHW)
hotspot<-ifelse(hotspot<1,hotspot==0,hotspot)

# Bind together the time vector and the hotspot vector
hotspot.20142016 <- cbind.data.frame(KBay.hotspot,hotspot)


###################################################################################
## Now calculate the Degree Heating Week (DHW) ##
# Set the window for the rolling sum (this is 24 half weeks = summing over 12 weeks)
k=24

# Create a dataframe including the time vector xi3
KBay.DHW <- data.frame(half.wk.dates)

# Calculate the DHW using rollapply, summing across a 12 week window (k=24 half weeks). Multiply value by 0.5 to get DHW
KBay.DHW$DHW<-(rollapply(hotspot.20142016$hotspot, width=k, by=1, sum, fill=c(NA,NA,NA), align="right",partial=FALSE)*0.5)

write.csv(KBay.DHW, "output/KBay.DHW.csv")

plot(KBay.DHW$half.wk.dates, KBay.DHW$DHW, type="l", xlab="Time", ylab=paste(expression("DHW (°C)")), col="red")
**Figure NA.** Degree Heating Weeks over 2014-2016 (data reported, figure not in mansucript).

Figure NA. Degree Heating Weeks over 2014-2016 (data reported, figure not in mansucript).

dev.copy(png, "output/KBay.DHWplot.png")
dev.off()

Temperature plots for logger temp at each site (color, top) and temperatures for the 2 bleachin years (2014-2015) (2015-2016). These plots, along with the benthic data (futher down) constitute Figure 1.

#######################
# Plot temp. daily mean temperatures for each reef
#######################
# Temperature and light data analysis ------
reefcolsBW <- c("gray60", "gray35")
reefcols2 <- c("mediumseagreen", "lightskyblue")

par(mfrow=c(1,1), mar=c(2,4,2,1), mgp=c(2,0.5,0))

#### Temperature plot
### NOAA-HIMB
k=1; lwd=0.7 # k-day moving averages
plot(mean ~ Date, annual.day.temps, type="l", ylab="Mean daily temp. (°C)", cex.lab=0.6, cex.axis=0.6, ylim=c(21, 31), xlab="", xaxt="n", col="gray85")
axis.Date(1, at=seq(min(annual.day.temps$Date), max(annual.day.temps$Date), by="3 mon"), 
          format="%b '%y", cex.axis=0.6, cex.lab=0.6)
legend("topright", lty=1, col=c("gray 80", reefcols2[c(1,2)]), legend=c("NOAA-HIMB", "Lilipuna", "Reef 14"), seg.len= 1.5, cex=0.6, lwd=1, bty="n", inset=c(0.01, 0), y.intersp = 1.3)

### Lilipuna
with(na.omit(data.frame(date=Lil_mean1$Date, mean=rollmean(Lil_mean1$mean, k, fill=NA))), {
  lines(date, mean, col=reefcols2[1], lwd=lwd)
})
with(na.omit(data.frame(date=Lil_mean2$Date, mean=rollmean(Lil_mean2$mean, k, fill=NA))), {
  lines(date, mean, col=reefcols2[1], lwd=lwd)
})
with(na.omit(data.frame(date=Lil_mean4$Date, mean=rollmean(Lil_mean4$mean, k, fill=NA))), { 
  lines(date, mean, col=reefcols2[1], lwd=lwd) 
})
# with(na.omit(data.frame(date=Lil_mean5$Date, mean=rollmean(Lil_mean5$mean, k, fill=NA))), { 
# lines(date, mean, col=reefcols2[1], lwd=1.2) })

### Reef 14
with(na.omit(data.frame(date=Rf14_mean1$Date, mean=rollmean(Rf14_mean1$mean, k, fill=NA))), { 
  lines(date, mean, col=reefcols2[2], lwd=lwd) 
})
with(na.omit(data.frame(date=Rf14_mean2$Date, mean=rollmean(Rf14_mean2$mean, k, fill=NA))), { 
  lines(date, mean, col=reefcols2[2], lwd=lwd) 
})
with(na.omit(data.frame(date=Rf14_mean3$Date, mean=rollmean(Rf14_mean3$mean, k, fill=NA))), { 
  lines(date, mean, col=reefcols2[2], lwd=lwd)
})
with(na.omit(data.frame(date=Rf14_mean4$Date, mean=rollmean(Rf14_mean4$mean, k, fill=NA))), { 
  lines(date, mean, col=reefcols2[2], lwd=lwd)
})
# with(na.omit(data.frame(date=Rf14_mean5$Date, mean=rollmean(Rf14_mean5$mean, k, fill=NA))), { 
# lines(date, mean, col=reefcols2[2], lwd=1.2) })

dev.copy(pdf, "figures/Temp.2014_2016_col.pdf", width=7, height=3)
dev.off()

temp.plot <- recordPlot() # records the plot for later


###########################
# Annual heating and cooling 
###########################

# Temperature and light data analysis ------

#### Temperature plot
# limit dates
Temp.event1<-annual.day.temps[(annual.day.temps$Date< "2015-03-01"),]
Temp.event2<-annual.day.temps[(annual.day.temps$Date>= "2015-01-01"),]

annual.color <- c("gray59", "black")
par(mfrow=c(1,1), mar=c(2,4,2,1), mgp=c(2,0.5,0))

### 2014-bleaching
k=1; lwd=0.7 # k-day moving averages
plot(mean ~ Date, Temp.event1, type="l", ylab="Mean daily temp. (°C)", cex.lab=0.6, cex.axis=0.6, ylim=c(21, 31), xlab="", xaxt="n", col=annual.color[1])
axis.Date(1, at=seq(min(Temp.event1$Date), max(Temp.event1$Date), by="3 mon"),
          format="%b", cex.axis=0.6, cex.lab=0.6)
abline(v=Temp.event1[Temp.event1$Date=="2014-10-10",1], col = "black", lwd=1, lty=2)
abline(v=Temp.event1[Temp.event1$Date=="2015-02-11",1], col = "black", lwd=1, lty=2)
legend("topleft", lty=1, col=annual.color, legend=c("2014-bleaching", "2015-bleaching"), seg.len= 1.5, cex=0.6, lwd=1, bty="n", inset=c(0.1, 0), y.intersp = 1.3)

### 2015 bleaching
par(new=T)
plot(mean ~ Date, Temp.event2, type="l", cex.lab=0.6, cex.axis=0.6, ylim=c(21, 31), xlab="", xaxt="n", ylab="", col=annual.color[2])
abline(v=Temp.event2[Temp.event2$Date=="2015-10-12",1], col = "gray80", lwd=1, lty=2) 
abline(v=Temp.event2[Temp.event2$Date=="2016-02-26",1], col = "gray80", lwd=1, lty=2)

dev.copy(pdf, "figures/Temp.AnnualHeating_col.pdf", width=7, height=3)
dev.off()
**Figure 1 (partial)**. Temperatures in Kāne'ohe Bay and the two reefs sites during repeated bleaching and recovery periods**Figure 1 (partial)**. Temperatures in Kāne'ohe Bay and the two reefs sites during repeated bleaching and recovery periods

Figure 1 (partial). Temperatures in Kāne’ohe Bay and the two reefs sites during repeated bleaching and recovery periods

Light

Irradiance at each reef was logged with loggers through time at ~1m depth. The data here was calibrated against a Li-Cor quantum sensor. As with temperature data, gaps represent logger failure or loss.

#setwd("data/environmental/Light data")
# Import light data

# Lilipuna
# 2014 Oct - 2016 Jun
Lil.PAR1<-rbind(
  read.csv("data/environmental/Light data/Lilipuna/SN_2485/Lilipuna_2485_T1_oct-dec 2014.csv"),
  read.csv("data/environmental/Light data/Lilipuna/SN_2485/Lilipuna_2485_T2_dec-feb 2015.csv"))

Lil.PAR2<-rbind(  
  read.csv("data/environmental/Light data/Lilipuna/SN_2485/Lilipuna_2485_T3_feb-apr 2015.csv"),
  read.csv("data/environmental/Light data/Lilipuna/SN_2485/Lilipuna_2485_T4_apr-jun 2015.csv"))

Lil.PAR4<-rbind(
  read.csv("data/environmental/Light data/Lilipuna/SN_2489/Lilipuna_2489_T6_sep-nov 2015.csv"),
  read.csv("data/environmental/Light data/Lilipuna/SN_2489/Lilipuna_2489_T7_nov-jan 2016.csv"),
  read.csv("data/environmental/Light data/Lilipuna/SN_2489/Lilipuna_2489_T8_jan-mar 2016.csv"))
  
# Reef 14
# 2014 Oct - 2016 Jun
Rf14.PAR1 <- rbind(
  read.csv("data/environmental/Light data/Reef 14/SN_2488/Reef 14_2488_T1_oct-dec 2014.csv"),
  read.csv("data/environmental/Light data/Reef 14/SN_2488/Reef 14_2488_T2_dec-feb 2015.csv"))

Rf14.PAR2 <- rbind(
  read.csv("data/environmental/Light data/Reef 14/SN_2488/Reef 14_2488_T3_feb-apr 2015.csv"),
  read.csv("data/environmental/Light data/Reef 14/SN_2488/Reef 14_2488_T4_apr-jun 2015.csv"))

Rf14.PAR3 <- 
  read.csv("data/environmental/Light data/Reef 14/SN_2488/Reef 14_2488_T5_jun-sep 2015.csv")

Rf14.PAR4 <- rbind(
  read.csv("data/environmental/Light data/Reef 14/SN_2488/Reef 14_2488_T6_sep-nov 2015.csv"),
  read.csv("data/environmental/Light data/Reef 14/SN_2488/Reef 14_2488_T7_nov-jan 2016.csv"),
  read.csv("data/environmental/Light data/Reef 14/SN_2488/Reef 14_2488_T8_jan-mar 2016.csv"))


# Set date to date class
Lil.PAR1$Date <- as.Date(Lil.PAR1$Date, format="%e/%m/%Y")
Lil.PAR2$Date <- as.Date(Lil.PAR2$Date, format="%e/%m/%Y")
Lil.PAR4$Date <- as.Date(Lil.PAR4$Date, format="%e/%m/%Y")

Rf14.PAR1$Date <- as.Date(Rf14.PAR1$Date, format="%e/%m/%Y")
Rf14.PAR2$Date <- as.Date(Rf14.PAR2$Date, format="%e/%m/%Y")
Rf14.PAR3$Date <- as.Date(Rf14.PAR3$Date, format="%e/%m/%Y")
Rf14.PAR4$Date <- as.Date(Rf14.PAR4$Date, format="%e/%m/%Y")

##########################
# Calibrate light data
# Lilipuna --- Logger SN: 2485 calibration: y = 23.674x - 48.339 
# Lilipuna --- Logger SN: 2489 calibration: y = 89.459x
# Reef 14  --- Logger SN: 2488 calibration: y = 74.242x
##########################

#Reef 14 -- Odyssey logger 2488

Rf14.PAR1<-Rf14.PAR1[,-5]; head(Rf14.PAR1) #remove "calibrated column"
# integrate raw values over time internal (15min * 60 sec) and calibrate to LiCor
Rf14.PAR1$LiCor_Calibrated<-(Rf14.PAR1$Raw*74.242/(15*60))
  
  Rf14.PAR2<-Rf14.PAR2[,-5]; head(Rf14.PAR2)
  Rf14.PAR2$LiCor_Calibrated<-(Rf14.PAR2$Raw*74.242/(15*60))
  
    Rf14.PAR3<-Rf14.PAR3[,-5]; head(Rf14.PAR3)
    Rf14.PAR3$LiCor_Calibrated<-(Rf14.PAR3$Raw*74.242/(15*60))
    
      Rf14.PAR4<-Rf14.PAR4[,-5]; head(Rf14.PAR4)
      Rf14.PAR4$LiCor_Calibrated<-(Rf14.PAR4$Raw*74.242/(15*60))
         

# Lilipuna -- Odyssey logger 2485 .....calibration
Lil.PAR1<-Lil.PAR1[,-5]; head(Lil.PAR1) 
Lil.PAR1$LiCor_Calibrated<-(Lil.PAR1$Raw*82.0/(15*60)); head(Lil.PAR1)
  Lil.PAR2<-Lil.PAR2[,-5]; head(Lil.PAR2)
  Lil.PAR2$LiCor_Calibrated<-(Lil.PAR2$Raw*82.0/(15*60)); head(Lil.PAR2)
  

# Lilipuna -- Odyssey logger 2489 calibration: 
Lil.PAR4<-Lil.PAR4[,-5]; head(Lil.PAR4)
Lil.PAR4$LiCor_Calibrated<-(Lil.PAR4$Raw*89.459/(15*60)); head(Lil.PAR4)

# Combine all calibrated data for Lilipuna, Reef 14
Lil.PAR  <- rbind(Lil.PAR1, Lil.PAR2, Lil.PAR4)

Rf14.PAR <- rbind(Rf14.PAR1, Rf14.PAR2, Rf14.PAR3, Rf14.PAR4)

#######################
#######################
# compiled data files

Lil.Temp  # Lilipuna temp data
Rf14.Temp # Reef 14 temp data

Lil.PAR   # Lilipuna light data
Rf14.PAR  # Reef 14 light data
########################

# Calculate daily light integrals
# Lilipuna
Lil.L1 <- aggregate(data.frame(mean=Lil.PAR1$LiCor_Calibrated), by=list(Date=Lil.PAR1$Date), FUN=mean)
Lil.L1$dli <- Lil.L1$mean * 0.0864 # Convert to mol.m2.d (daily light integral)

Lil.L2 <- aggregate(data.frame(mean=Lil.PAR2$LiCor_Calibrated), by=list(Date=Lil.PAR2$Date), FUN=mean); Lil.L2$dli <- Lil.L2$mean * 0.0864

Lil.L4 <- aggregate(data.frame(mean=Lil.PAR4$LiCor_Calibrated), by=list(Date=Lil.PAR4$Date), FUN=mean); Lil.L4$dli <- Lil.L4$mean * 0.0864
Lil.L4<-Lil.L4[!(Lil.L4$Date=="2016-03-01"),]


#Reef 14
Rf14.L1 <- aggregate(data.frame(mean=Rf14.PAR1$LiCor_Calibrated), by=list(Date=Rf14.PAR1$Date), FUN=mean); Rf14.L1$dli <- Rf14.L1$mean * 0.0864

Rf14.L2 <- aggregate(data.frame(mean=Rf14.PAR2$LiCor_Calibrated), by=list(Date=Rf14.PAR2$Date), FUN=mean); Rf14.L2$dli <- Rf14.L2$mean * 0.0864

Rf14.L3 <- aggregate(data.frame(mean=Rf14.PAR3$LiCor_Calibrated), by=list(Date=Rf14.PAR3$Date), FUN=mean); Rf14.L3$dli <- Rf14.L3$mean * 0.0864

Rf14.L4 <- aggregate(data.frame(mean=Rf14.PAR4$LiCor_Calibrated), by=list(Date=Rf14.PAR4$Date), FUN=mean); Rf14.L4$dli <- Rf14.L4$mean * 0.0864


# testing for differences among reefs

## DAILY INTEGRATED LIGHT
# Light data dli
# Lil.PAR1-2,4-5  and Rf14.PAR1-5

Lil.PAR.dli.all<-rbind(Lil.L1, Lil.L2, Lil.L4)
Rf14.PAR.dli.all<-rbind(Rf14.L1, Rf14.L2, Rf14.L3, Rf14.L4)



###########################
# Plot daily light integrals
###########################

# Light data  ------
reefcolsBW <- c("gray60", "gray35")
reefcols2 <- c("mediumseagreen", "lightskyblue")

par(mfrow=c(1,1), mar=c(2,4,2,1), mgp=c(2,0.5,0))
k=1; lwd=0.8 # k-day moving averages
plot(mean ~ Date, Lil.PAR.dli.all, type="n", ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), xaxt="n", xlab="", cex.lab=0.7, ylim=c(0,50), cex.axis=0.7)
axis.Date(1, at=seq(min(Lil.PAR.dli.all$Date), max(Lil.PAR.dli.all$Date), by="1 mon"), format="%b '%y", cex.axis=0.7, cex.lab=0.7)
with(na.omit(data.frame(date=Lil.L1$Date, dli=rollmean(Lil.L1$dli, k, fill=NA))), lines(date, dli, col=reefcols2[1], lwd=lwd))
with(na.omit(data.frame(date=Lil.L2$Date, dli=rollmean(Lil.L2$dli, k, fill=NA))), lines(date, dli, col=reefcols2[1], lwd=lwd))
with(na.omit(data.frame(date=Lil.L4$Date, dli=rollmean(Lil.L4$dli, k, fill=NA))), lines(date, dli, col=reefcols2[1], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L1$Date, dli=rollmean(Rf14.L1$dli, k, fill=NA))), lines(date, dli, col=reefcols2[2], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L2$Date, dli=rollmean(Rf14.L2$dli, k, fill=NA))), lines(date, dli, col=reefcols2[2], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L3$Date, dli=rollmean(Rf14.L3$dli, k, fill=NA))), lines(date, dli, col=reefcols2[2], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L4$Date, dli=rollmean(Rf14.L4$dli, k, fill=NA))), lines(date, dli, col=reefcols2[2], lwd=lwd))
legend("topright", lty=1, col=c(reefcols2[c(1,2)]), legend=c("Lilipuna", "Reef 14"), seg.len= 1.2, cex=0.8, y.intersp=1.2, x.intersp = 0.8, lwd=2, bty="n")
**Supplementary Figure 1.**  Light conditions at the two reefs through time.

Supplementary Figure 1. Light conditions at the two reefs through time.

dev.copy(pdf, "figures/PAR.pdf", width=7, height=3)
dev.off()

Ecological data

2014-2016 benthic

#### Bleaching and recovery for year 1 and 2
ecology<-read.csv("data/ecology/Reef survey 2014-2016.csv") #import your data file here

# remove "Slope" and Flat habitat where corals were not collected
ecology<-ecology[(ecology$Habitat=="Crest") ,] 
droplevels(ecology$Habitat)

ecology$total.MC.bleach<-ecology$MC.pale+ecology$MC.white

##### make a dataframe for data means by transect (= unit of replication)

coral<-aggregate(coral.cover~Transect+Site+Status+Period, data=ecology, sum)
healthy<-aggregate(coral.pigmented~Transect+Site+Status+Period, data=ecology, sum)
pale<-aggregate(coral.pale~Transect+Site+Status+Period, data=ecology, sum)
white<-aggregate(coral.white~Transect+Site+Status+Period, data=ecology, sum)
MC<-aggregate(MC~Transect+Site+Status+Period, data=ecology, sum)
MC.healthy<-aggregate(MC.pigmented~Transect+Site+Status+Period, data=ecology, sum)
MC.pale<-aggregate(MC.pale~Transect+Site+Status+Period, data=ecology, sum)
MC.white<-aggregate(MC.white~Transect+Site+Status+Period, data=ecology, sum)

df.ecol<-cbind(coral, healthy[c(5,0)], pale[c(5,0)], white[c(5,0)], MC[c(5,0)], MC.healthy[c(5,0)], MC.pale[c(5,0)], MC.white[c(5,0)])
df.ecol<- df.ecol %>% mutate(
  totalcoral.bleach= coral.pale + coral.white,
  totalMC.bleach= MC.pale + MC.white)

eco.surv<-subset(df.ecol, select = c(-coral.pale, -coral.white, -MC.pale, -MC.white))
  
# use "eco.surv" to test mean and SD 
# df.m + dfSE to make figures

coral<-aggregate(coral.cover~Site+Status+Period, data=eco.surv, mean)
healthy<-aggregate(coral.pigmented~Site+Status+Period, data=eco.surv, mean)
bleached<-aggregate(totalcoral.bleach~Site+Status+Period, data=eco.surv, mean)

MC<-aggregate(MC~Site+Status+Period, data=eco.surv, mean)
MC.healthy<-aggregate(MC.pigmented~Site+Status+Period, data=eco.surv, mean)
MC.bleach<-aggregate(totalMC.bleach~Site+Status+Period, data=eco.surv, mean)

df.ecol.m<-cbind(coral, healthy[c(4,0)], bleached[c(4,0)], MC[c(4,0)], MC.healthy[c(4,0)], MC.bleach[c(4,0)])

# SD
coralSD<-aggregate(coral.cover~Site+Status+Period, data=eco.surv, sd)
healthySD<-aggregate(coral.pigmented~Site+Status+Period, data=eco.surv, sd)
bleachedSD<-aggregate(totalcoral.bleach~Site+Status+Period, data=eco.surv, sd)

MCSD<-aggregate(MC~Site+Status+Period, data=eco.surv, sd)
MC.healthySD<-aggregate(MC.pigmented~Site+Status+Period, data=eco.surv, sd)
MC.bleachSD<-aggregate(totalMC.bleach~Site+Status+Period, data=eco.surv, sd)

df.ecol.SD<-cbind(coralSD, healthySD[c(4,0)], bleachedSD[c(4,0)], MCSD[c(4,0)], MC.healthySD[c(4,0)], MC.bleachSD[c(4,0)])


df.fig<-cbind(df.ecol.m, df.ecol.SD[c(4:9,0)]); colnames(df.fig) <-c("Site", "Status", "Period", "coral.cover", "coral.pigmented", "totalcoral.bleach", "MC", "MC.pigmented", "totalMC.bleach", "coral.SD", "coral.pigSD", "coral.blSD", "MCSD", "MC.pigSD", "MC.blSD"); df.fig

df.fig<-df.fig %>% mutate(
  groupyear = ifelse(Period =="2014 Oct" | Period =="2015 Feb",
                     "2014-2015", "2015-2016"))

########################
## color palettes ######
########################
reefcols <- c("gray60", "gray35")
ecol.labs<-c("Oct'14 \nBl","Feb'15 \nRc","Oct'15 \nBl", "Feb'16 \nRc")

text.formatting<-theme_classic() + 
  theme(text=element_text(size=8),
  axis.text.y=element_text(
      margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=8), 
    axis.text.x=element_text(
      margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=8), 
  plot.title = element_text(hjust = 0.5))

pd <- position_dodge(0.15) #offset for error bars

##############
# MC corals
FigMC <- ggplot(df.fig, aes(x=Period, y=MC, group=Site, color=Site)) + 
  geom_errorbar(aes(ymin=MC-MCSD, ymax= MC+MCSD), width=0, position=pd) +
  geom_line(position=pd, size=.6) +
  geom_point(position=pd, size=1) +
  ylab("") +
  scale_color_manual(values=reefcols2) +
  ggtitle(expression(bold(paste(~bolditalic("M. capitata")~cover, sep="")))) +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(0,1.05),oob = rescale_none)+ 
  scale_x_discrete(name ="Year and Period", 
      labels=ecol.labs)+
  text.formatting

# bleached MC coral
FigBlMC <- ggplot(df.fig, aes(x=Period, y=totalMC.bleach, group=Site, color=Site)) + 
  geom_errorbar(aes(ymin=totalMC.bleach-MC.blSD, ymax= totalMC.bleach+MC.blSD), width=0.0, position=pd) +
  geom_line(position=pd, size=.6) +
  geom_point(position=pd, size=1) +
  ylab("") +
  xlab(ecol.labs) +
  scale_color_manual(values=reefcols2) +
  ggtitle(expression(bold(paste(~Bleached, bolditalic(" M. capitata")~cover, sep="")))) +
  scale_y_continuous(limits=c(0,1.05),oob = rescale_none)+ 
  scale_x_discrete(name ="Year and Period", 
      labels=ecol.labs)+ text.formatting

############
############
#export figures

#############
# color figures all corals
Ecol.Fig<-plot_grid(FigMC, FigBlMC,
     labels=c('b', 'c'), label_size=8, hjust=-3, vjust= 3, ncol=2, nrow=1);
Ecol.Fig
**Figure 2 (partial)**. *Montipora capitata* total cover and bleached cover in Kāne'ohe Bay at the two reefs sites during repeated bleaching and recovery periods.

Figure 2 (partial). Montipora capitata total cover and bleached cover in Kāne’ohe Bay at the two reefs sites during repeated bleaching and recovery periods.

##### save the figure and export to directory? ####
dev.copy(pdf, "figures/Ecol.Fig_col.pdf", width=7, height=4)
dev.off()

Statistics for bleaching and recovery bethic data.

The first table output shows the effects of bleaching on total Montipora cover and post-hoc test for the 4 periods, ‘sliced’ by Period to show Site effects.

##################
###############
##########
#####
#
#  CORAL COMMUNITY BLEACHING ACROSS SITES

# dataframe is 'ecology'

################
# total MC cover
################

MC.cover.mod<-lm(asin(sqrt(MC))~Period*Site, data=ecology)
anova(MC.cover.mod)
## Analysis of Variance Table
## 
## Response: asin(sqrt(MC))
##              Df Sum Sq  Mean Sq F value    Pr(>F)    
## Period        3 0.2954 0.098461  7.4871 7.445e-05 ***
## Site          1 0.0248 0.024838  1.8887 0.1703381    
## Period:Site   3 0.2803 0.093426  7.1042 0.0001245 ***
## Residuals   312 4.1031 0.013151                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(MC.cover.mod, ~Site|Period)
contrast(posthoc)
## Period = 2014 Oct:
##  contrast        estimate     SE  df t.ratio p.value
##  Lilipuna effect -0.05582 0.0128 312 -4.354  <.0001 
##  Reef 14 effect   0.05582 0.0128 312  4.354  <.0001 
## 
## Period = 2015 Feb:
##  contrast        estimate     SE  df t.ratio p.value
##  Lilipuna effect  0.02308 0.0128 312  1.801  0.0727 
##  Reef 14 effect  -0.02308 0.0128 312 -1.801  0.0727 
## 
## Period = 2015 Oct:
##  contrast        estimate     SE  df t.ratio p.value
##  Lilipuna effect  0.00774 0.0128 312  0.604  0.5465 
##  Reef 14 effect  -0.00774 0.0128 312 -0.604  0.5465 
## 
## Period = 2016 Feb:
##  contrast        estimate     SE  df t.ratio p.value
##  Lilipuna effect -0.01024 0.0128 312 -0.799  0.4250 
##  Reef 14 effect   0.01024 0.0128 312  0.799  0.4250 
## 
## Note: contrasts are still on the asin.sqrt scale 
## P value adjustment: fdr method for 2 tests
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna 0.0715 0.0181 312  0.03587   0.1072  a    
##  Reef 14  0.1832 0.0181 312  0.14752   0.2189   b   
## 
## Period = 2015 Feb:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Reef 14  0.0271 0.0181 312 -0.00862   0.0627  a    
##  Lilipuna 0.0732 0.0181 312  0.03755   0.1089  a    
## 
## Period = 2015 Oct:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Reef 14  0.0499 0.0181 312  0.01418   0.0855  a    
##  Lilipuna 0.0653 0.0181 312  0.02966   0.1010  a    
## 
## Period = 2016 Feb:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna 0.0773 0.0181 312  0.04161   0.1130  a    
##  Reef 14  0.0978 0.0181 312  0.06209   0.1334  a    
## 
## Results are given on the asin(sqrt(mu)) (not the response) scale. 
## Confidence level used: 0.95 
## Note: contrasts are still on the asin.sqrt scale 
## significance level used: alpha = 0.05

The second table output shows the effects of bleaching on BLEACHED Montipora cover and post-hoc test for the 4 periods. THe post-hoc tests show effects of Site, then the effects of Site within each period.

################
# bleached MC
################

MC.bleach.mod<-lm(total.MC.bleach~Period*Site, data=ecology)
anova(MC.bleach.mod)
## Analysis of Variance Table
## 
## Response: total.MC.bleach
##              Df  Sum Sq   Mean Sq F value   Pr(>F)   
## Period        3 0.01540 0.0051331  2.4232 0.065855 . 
## Site          1 0.02059 0.0205868  9.7186 0.001994 **
## Period:Site   3 0.00215 0.0007164  0.3382 0.797716   
## Residuals   312 0.66091 0.0021183                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(MC.bleach.mod, ~Site)
contrast(posthoc)
##  contrast        estimate      SE  df t.ratio p.value
##  Lilipuna effect -0.00802 0.00257 312 -3.117  0.0020 
##  Reef 14 effect   0.00802 0.00257 312  3.117  0.0020 
## 
## Results are averaged over the levels of: Period 
## P value adjustment: fdr method for 2 tests
multcomp::cld(posthoc, Letters=letters)
##  Site      lsmean      SE  df lower.CL upper.CL .group
##  Lilipuna 0.00438 0.00364 312 -0.00278   0.0115  a    
##  Reef 14  0.02042 0.00364 312  0.01326   0.0276   b   
## 
## Results are averaged over the levels of: Period 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(MC.bleach.mod, ~Site|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site      lsmean      SE  df  lower.CL upper.CL .group
##  Lilipuna 0.01250 0.00728 312 -0.001819   0.0268  a    
##  Reef 14  0.03083 0.00728 312  0.016515   0.0452  a    
## 
## Period = 2015 Feb:
##  Site      lsmean      SE  df  lower.CL upper.CL .group
##  Lilipuna 0.00000 0.00728 312 -0.014319   0.0143  a    
##  Reef 14  0.00833 0.00728 312 -0.005985   0.0227  a    
## 
## Period = 2015 Oct:
##  Site      lsmean      SE  df  lower.CL upper.CL .group
##  Lilipuna 0.00500 0.00728 312 -0.009319   0.0193  a    
##  Reef 14  0.02750 0.00728 312  0.013181   0.0418   b   
## 
## Period = 2016 Feb:
##  Site      lsmean      SE  df  lower.CL upper.CL .group
##  Lilipuna 0.00000 0.00728 312 -0.014319   0.0143  a    
##  Reef 14  0.01500 0.00728 312  0.000681   0.0293  a    
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05

Physiology-Immunity

Import data and normalize
Raw data from physiological assessments is imported and normalized according to established protocols. Immunolgical data has been normalized previously and is represented as the final dataframe.

###########################################################################
## data period: October 2014, February 2015, October 2015, February 2016 ##
###########################################################################

### data file
# all lab data from all time points (physio and immuno, PLUS color analysis)
physimmun<-read.csv("data/Gates_Mydlarz_20142016_physimmun.csv", header=T, na.string=NA) #physimmuno data
qPCR.data<-read.csv("data/Gates_Mydlarz_20142016_qPCR.csv", header=T, na.string=NA) #qPCR data

data<-read.csv("data/Gates_Mydlarz_20142016_ALL_DATA.csv", header=T, na.string=NA) #qPCR data


################################################
# PHYSIOLOGY: calculate, normalize dependent variables
################################################
data$cells..ml<-as.numeric(data$cells..ml)
data$ID<-as.factor(data$ID)

# helpful shorthand
SA<-data$surface.area.cm2 # surface area in cm2
blastate<-data$total.blastate.ml # tissue slurry blastate in ml

# Symbiodinium.cells..cm2 == cell.ml * blastate / SA
data$zoox<- (data$cells..ml*blastate)/SA

# ug.chlorophyll.a..cm2 == ug.chl.a.ml * blastate / SA
data$chla<- (data$ug.chla..ml*blastate)/SA

# pg.chlorophyll.a..cell == ug.chl.a.ml * 10^6 / cells.ml
data$chlacell<- (data$ug.chla..ml*10^6)/data$cells..ml

# AFDW.mg..cm2 == convert AFDW g to mg, mutiply by blastate volume, divide by cm2
data$AFDW<- (data$g.AFDW..ml*1000*blastate)/SA

# mg.protein..cm2 == mg.protein.ml * blastate / SA
data$prot<- (data$mg.prot..ml*blastate)/SA


#### rearrange and clean up

# drop unwanted columns by name
data<-data[!colnames(data) %in% c("total.blastate.ml", "surface.area.cm2", "mg.prot..ml", "cells..ml", "ug.chla..ml", "g.AFDW..ml")]

data<-data[!(data$ID=="N37"), ] # this coral has an NA for symbiont type, remove from df

# reorder columns to show: hierarchy of structure, physio., immuno., colorimetry
# this is now the working dataframe
df<-data[, c(1:6, 16:18,20,19, 7:15)]

# remove negative values for MEL and POX and
df$MEL[df$MEL <= 0]<- NA
df$MEL[df$MEL > 0.05]<- NA

df$POX[df$POX <= 0]<- NA
df$CAT[df$CAT <= 0]<- NA
df$PPO[df$PPO <= 0]<- NA


df$chlacell[df$chlacell >= 10]<- NA # removing 2 outliers
df$CAT[df$CAT >= 900]<- NA # 2 outliers removed

# don't need "Ramp" lab data, this was lab thermally stressed. Remove this.
df<-df[!(df$Status=="Lab-Ramp"),]


#########################################################
#########################################################
#########################################################

### Summary dataframes
# These dataframes show mean, SE and n for each group. This can be subset to make figures or exported as a summary file. The header of the dataframe is shown below...

# make summary dataframes

# dataframes of means can later be divided into categories by "C or D" dominated
###------#### means

# physiology
zoox.m<-aggregate(zoox~Site+Status+Period+dom,data=df, mean)
chla.m<-aggregate(chla~Site+Status+Period+dom,data=df, mean)
chlacell.m<-aggregate(chlacell~Site+Status+Period+dom, data=df, mean)
prot.m<-aggregate(prot~Site+Status+Period+dom, data=df, mean)
AFDW.m<-aggregate(AFDW~Site+Status+Period+dom, data=df, mean)

# imunology
CAT.m<-aggregate(CAT~Site+Status+Period+dom,data=df, mean)
POX.m<-aggregate(POX~Site+Status+Period+dom,data=df, mean)
SOD.m<-aggregate(SOD~Site+Status+Period+dom,data=df, mean)
PPO.m<-aggregate(PPO~Site+Status+Period+dom,data=df, mean)
MEL.m<-aggregate(MEL~Site+Status+Period+dom,data=df, mean)

###------#### SE
# physiology
zoox.SE<-aggregate(zoox~Site+Status+Period+dom,data=df, std.error)
chla.SE<-aggregate(chla~Site+Status+Period+dom,data=df, std.error)
chlacell.SE<-aggregate(chlacell~Site+Status+Period+dom, data=df, std.error)
prot.SE<-aggregate(prot~Site+Status+Period+dom, data=df, std.error)
AFDW.SE<-aggregate(AFDW~Site+Status+Period+dom, data=df, std.error)

# imunology
CAT.SE<-aggregate(CAT~Site+Status+Period+dom,data=df, std.error)
POX.SE<-aggregate(POX~Site+Status+Period+dom,data=df, std.error)
SOD.SE<-aggregate(SOD~Site+Status+Period+dom,data=df, std.error)
PPO.SE<-aggregate(PPO~Site+Status+Period+dom,data=df, std.error)
MEL.SE<-aggregate(MEL~Site+Status+Period+dom,data=df, std.error)

###------#### n-value
# physiology
zoox.n<-aggregate(zoox~Site+Status+Period+dom,data=df, length)
chla.n<-aggregate(chla~Site+Status+Period+dom,data=df, length)
chlacell.n<-aggregate(chlacell~Site+Status+Period+dom, data=df, length)
prot.n<-aggregate(prot~Site+Status+Period+dom, data=df, length)
AFDW.n<-aggregate(AFDW~Site+Status+Period+dom, data=df, length)

# imunology
CAT.n<-aggregate(CAT~Site+Status+Period+dom,data=df, length)
POX.n<-aggregate(POX~Site+Status+Period+dom,data=df, length)
SOD.n<-aggregate(SOD~Site+Status+Period+dom,data=df, length)
PPO.n<-aggregate(PPO~Site+Status+Period+dom,data=df, length)
MEL.n<-aggregate(MEL~Site+Status+Period+dom,data=df, length)


#### physio dataframe
fig.df.phys<-cbind(zoox.m, zoox.SE[c(5,0)],chla.m[c(5,0)], chla.SE[c(5,0)], chlacell.m[c(5,0)], chlacell.SE[c(5,0)], prot.m[c(5,0)], prot.SE[c(5,0)], AFDW.m[c(5,0)], AFDW.SE[c(5,0)], chla.n[c(5,0)])

colnames(fig.df.phys)<-c("Site","Status", "Period", "dom", "zoox", "zoox.SE", "chla", "chla.SE", "chlacell", "chlacell.SE", "prot", "prot.SE", "AFDW", "AFDW.SE", "N-chla")

# make an interaction column
fig.df.phys$int<-interaction(fig.df.phys$dom, fig.df.phys$Site)
fig.df.phys$int<-factor(fig.df.phys$int)
fig.df.phys$SimpleSite<-c("Lilipuna", "Reef 14")
####
#### 

#### immuno dataframe
fig.df.immuno<-cbind(CAT.m, CAT.SE[c(5,0)], POX.m[c(5,0)], POX.SE[c(5,0)], SOD.m[c(5,0)], SOD.SE[c(5,0)], PPO.m[c(5,0)], PPO.SE[c(5,0)], MEL.m[c(5,0)], MEL.SE[c(5,0)])

colnames(fig.df.immuno)<-c("Site","Status", "Period", "dom", "CAT", "CAT.SE", "POX", "POX.SE", "SOD", "SOD.SE", "PPO", "PPO.SE", "MEL", "MEL.SE")


#### Use this code for 'PRE' Site immunity data without symbiont community data
CAT.m2<-aggregate(CAT~Site+Status+Period,data=df, mean)
POX.m2<-aggregate(POX~Site+Status+Period,data=df, mean)
SOD.m2<-aggregate(SOD~Site+Status+Period,data=df, mean)
PPO.m2<-aggregate(PPO~Site+Status+Period,data=df, mean)
MEL.m2<-aggregate(MEL~Site+Status+Period,data=df, mean)

### PRE------#### SE immunity
CAT.SE2<-aggregate(CAT~Site+Status+Period,data=df, std.error)
POX.SE2<-aggregate(POX~Site+Status+Period,data=df, std.error)
SOD.SE2<-aggregate(SOD~Site+Status+Period,data=df, std.error)
PPO.SE2<-aggregate(PPO~Site+Status+Period,data=df, std.error)
MEL.SE2<-aggregate(MEL~Site+Status+Period,data=df, std.error)

### PRE------#### n immunity (n = 6-8 for PRE-field)
CAT.2n<-aggregate(CAT~Site+Status+Period,data=df, length)
POX.2n<-aggregate(POX~Site+Status+Period,data=df, length)
SOD.2n<-aggregate(SOD~Site+Status+Period,data=df, length)
PPO.2n<-aggregate(PPO~Site+Status+Period,data=df, length)
MEL.2n<-aggregate(MEL~Site+Status+Period,data=df, length)

Pre.immuno<-cbind(CAT.m2[c(1:2),], CAT.SE2[c(1:2),4], POX.m2[c(1:2),4], POX.SE2[c(1:2),4], SOD.m2[c(1:2),4], SOD.SE2[c(1:2),4], PPO.m2[c(1:2),4], PPO.SE2[c(1:2),4], MEL.m2[c(1:2),4], MEL.SE2[c(1:2),4])
colnames(Pre.immuno)<-c("Site","Status", "Period", "CAT", "CAT.SE", "POX", "POX.SE", "SOD", "SOD.SE", "PPO", "PPO.SE", "MEL", "MEL.SE")

##########
Pre.immuno$dom<-"unident"
Pre.immuno<-Pre.immuno[, c(1:3, 14, 4:13)]

## now this subset of site means can be overlayed with datafame from Oct 2014- Feb 2016
fig.df.immuno.comp<-rbind(Pre.immuno, fig.df.immuno)

# make an interaction column
fig.df.immuno.comp$int<-interaction(fig.df.immuno.comp$dom, fig.df.immuno.comp$Site)
fig.df.immuno.comp$int<-factor(fig.df.immuno.comp$int)
fig.df.immuno.comp$SimpleSite<-c("Lilipuna", "Reef 14")

fig.df.immuno.comp$int<-factor(fig.df.immuno.comp$int, levels=c("unident.Lilipuna.Pre", "unident.Reef 14.Pre", "C.Lilipuna", "C.Reef 14", "D.Lilipuna", "D.Reef 14"))

Multivariate analyses

Multivariate non-dimensional scaling (NMDS) plots for the coral physiotype trajectories and the plot of physiotypes through time with response vectors.

# multivariate test using adonis-- Manova, mutivariate analysis of variance
df.MV<-df

# remove columns unnecessary for final analysis, few factors retained
df.MV<-df.MV[ , !names(df.MV) %in% c("Species", "Status_Site", "ID", "propC", "propD", "syms")]

df.MV$Period.dom<-interaction(df.MV$Period, df.MV$dom)

# remove prior-to-bleaching timepoints and unused levels that arise from this action
df.MV<-df.MV[(!df.MV$Period=="2014 Lab"),] 
df.MV<-df.MV[(!df.MV$Period=="2014 Field.Feb"),]
df.MV$Period<-droplevels(df.MV$Period) 
df.MV$Status<-droplevels(df.MV$Status)
df.MV$Period.dom<-droplevels(df.MV$Period.dom)

#drop all NAs, removing rows where NA was observed
df.MV.noNA<-df.MV %>% drop_na()

Coral trajectories

The trajectories of corals in multivariate traitspace. Here is the single panel NMDS for Lilipuna and Reef 14 corals dominated by C- or D-symbionts through the 4 periods of bleaching and recovery.

PERMANOVA output here is from adonis2 with a scaled and cnetered matrix.

# cleaned data for full dataframe NMDS (Site, Period, dom)

fullNMDS<-df.MV.noNA
fullNMDS$Period.Site.dom<-interaction(fullNMDS$Period, fullNMDS$Site, fullNMDS$dom)

fac.fullNMDS<-fullNMDS %>% 
  select(c(Period, Status, Site, dom, Period.Site.dom)) #  sans responses

resp.fullNMDS<-fullNMDS[, 4:13] # just responses

resp.fullNMDS.scaled <- scale(resp.fullNMDS, center = TRUE, scale = TRUE) 

##  Run MANOVA
all.Manova<-adonis2(resp.fullNMDS.scaled~Period*Site*dom, data=fullNMDS, permutations=1000,  method="euclidian")
all.Manova
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = resp.fullNMDS.scaled ~ Period * Site * dom, data = fullNMDS, permutations = 1000, method = "euclidian")
##                  Df SumOfSqs      R2       F   Pr(>F)    
## Period            3   881.55 0.31041 50.8698 0.000999 ***
## Site              1    61.45 0.02164 10.6376 0.000999 ***
## dom               1   196.12 0.06906 33.9511 0.000999 ***
## Period:Site       3    65.49 0.02306  3.7790 0.000999 ***
## Period:dom        3    55.44 0.01952  3.1989 0.000999 ***
## Site:dom          1     5.25 0.00185  0.9089 0.485514    
## Period:Site:dom   3    20.82 0.00733  1.2015 0.229770    
## Residual        269  1553.88 0.54714                     
## Total           284  2840.00 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## **not used**
##### test dispersion using PERMDISP
disp <- betadisper(vegdist(resp.fullNMDS.scaled, method = "euclidian"), fac.fullNMDS$Period)

#test for overall differences
#anova(disp)

## Permutation test for pairwise comparisons
#permutest(disp, pairwise = TRUE)

NMDS trait space trajectory First showing the observed dissimilarily and fit of NMDS in stress lot

###
allNMDS<-metaMDS(resp.fullNMDS, distance='bray', k=2, trymax=100) 
stressplot(allNMDS)

…now the 4 part convex hulls for coral trajectories.

allNMDS.df<-data.frame(x=allNMDS$point[,1], y=allNMDS$point[,2], 
                            Period=as.factor(fac.fullNMDS[,1]),
                            Site=as.factor(fac.fullNMDS[,2]),
                            Status=as.factor(fac.fullNMDS[,3]),
                            dom=as.factor(fac.fullNMDS[,4]),
                            Period.Site.dom=as.factor(fac.fullNMDS[,5]))

colnames(allNMDS.df)[1:2]<-c("MDS1", "MDS2")

#centroid means
NMDS.mean=aggregate(allNMDS.df[,1:2],list(group=allNMDS.df$Period.Site.dom), mean)

############### All groups NMDS
fullNMDS$Period.Site.dom<-droplevels(fullNMDS$Period.Site.dom)
groups<-fullNMDS$Period.Site.dom
groups2<-c("C-Lilipuna", "C-Reef 14", "D-Lilipuna", "D-Reef 14")
colors=c(
  "firebrick1", #C Lilipuna
  "firebrick3", #C Reef 14
  "steelblue1", # D Lilipuna
  "steelblue3") # D Reef 14

### 

plot<-ordiplot(allNMDS, type="n", cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.3, 0.3))
axis(side = 1, labels = FALSE)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
ordihull(plot, groups=fullNMDS$Period, draw="polygon", alpha=30, col=c("gray20", "gray40", "gray30", "gray50"), border=FALSE)

# points for start and end
with(plot, points(NMDS.mean[1,2], NMDS.mean[1,3], col=colors[1], pch=16, cex=0.7, lwd=1)) # Lil C
with(plot, points(NMDS.mean[4,2], NMDS.mean[4,3], col=colors[1], pch=17, cex=1, lwd=1)) # Lil C
with(plot, points(NMDS.mean[5,2], NMDS.mean[5,3], col=colors[2], pch=16, cex=0.7, lwd=1)) # R14 C
with(plot, points(NMDS.mean[8,2], NMDS.mean[8,3], col=colors[2], pch=17, cex=1, lwd=1)) # R14 C
with(plot, points(NMDS.mean[9,2], NMDS.mean[9,3], col=colors[3], pch=16, cex=0.7, lwd=1)) # Lil D
with(plot, points(NMDS.mean[12,2], NMDS.mean[12,3], col=colors[3], pch=17, cex=1, lwd=1)) # Lil D
with(plot, points(NMDS.mean[13,2], NMDS.mean[13,3], col=colors[4], pch=16, cex=0.7, lwd=1)) # R14 D
with(plot, points(NMDS.mean[16,2], NMDS.mean[16,3], col=colors[4], pch=17, cex=1, lwd=1)) # R14 D

# points for others (if want)
with(plot, points(NMDS.mean[2,2], NMDS.mean[2,3], col=colors[1], pch=16, cex=0.7, lwd=1)) # Lil C
with(plot, points(NMDS.mean[3,2], NMDS.mean[3,3], col=colors[1], pch=16, cex=0.7, lwd=1)) # Lil C
with(plot, points(NMDS.mean[6,2], NMDS.mean[6,3], col=colors[2], pch=16, cex=0.7, lwd=1)) # R14 C
with(plot, points(NMDS.mean[7,2], NMDS.mean[7,3], col=colors[2], pch=16, cex=0.7, lwd=1)) # R14 C
with(plot, points(NMDS.mean[10,2], NMDS.mean[10,3], col=colors[3], pch=16, cex=0.7, lwd=1)) # Lil D
with(plot, points(NMDS.mean[11,2], NMDS.mean[11,3], col=colors[3], pch=16, cex=0.7, lwd=1)) # Lil D
with(plot, points(NMDS.mean[14,2], NMDS.mean[14,3], col=colors[4], pch=16, cex=0.7, lwd=1)) # R14 D
with(plot, points(NMDS.mean[15,2], NMDS.mean[15,3], col=colors[4], pch=16, cex=0.7, lwd=1)) # R14 D

# lines for start and end
with(plot, lines(NMDS.mean[c(1:4),2], NMDS.mean[c(1:4),3], lwd=1.5, col =colors[1], lty=2))# Lil C
with(plot, lines(NMDS.mean[c(5:8),2], NMDS.mean[c(5:8),3], lwd=1.5, col =colors[2], lty=1))# R14 C
with(plot, lines(NMDS.mean[c(9:12),2], NMDS.mean[c(9:12),3], lwd=1.5, col =colors[3], lty=2))# Lil D
with(plot, lines(NMDS.mean[c(13:16),2], NMDS.mean[c(13:16),3], lwd=1.5, col =colors[4], lty=1))# R14 D

# add labels
text(-0.38, 0.28, "Bleaching \n 2015", cex=0.7, col="gray50")
text(-0.28, -0.2, "Bleaching \n 2014", cex=0.7, col="gray50")
text(0.25, -0.1, "Recovery \n 2015", cex=0.7, col="gray50")
text(0.07, 0.21, "Recovery \n 2016", cex=0.7, col="gray50")
legend("topright", legend=groups2, inset=c(0.03, -0.01), x.intersp=0.6, y.intersp=0.9, cex=0.7, pch=16, col=colors, lty=c(1,3,1,3), seg.len=2, pt.cex=1, bty="n")
**Figure 3.** Non-metric multidimensional scaling (NMDS) analyses of coral performance envelopes during bleaching stress and post-bleaching recovery.

Figure 3. Non-metric multidimensional scaling (NMDS) analyses of coral performance envelopes during bleaching stress and post-bleaching recovery.

dev.copy(pdf, "figures/1panel_allNMDS.pdf", height=6, width=6, useDingbats=FALSE)
dev.off() 

Coral physiotypes

The physiotypes of corals in multivariate traitspace. Here are the four panel NMDS for Lilipuna and Reef 14 Physiology and Immunity data in first bleaching-recovery (2014-2015) and second bleaching-recovery (2015-2016).

code for NMDS analysis separated by sites in each time period.

########## ########## ########## ########## ########## 
 ########## ########## ########## ########## ########## 
########## ########## ########## ########## ########## 

# for 4 panel NMDS separated by site and B-R periods
## build Manova dataframes 

########## Lilipuna 2014-2015
## Lilipuna only #######

## all data with lilupuna only
Manova.Lil.lev.dat<-df.MV.noNA %>% # "lev.dat" has levels and data
  filter(Site=="Lilipuna") 

# lilipuna bleaching-recovery1
Manova.Lil.BR1.lev.dat<-Manova.Lil.lev.dat %>% # "BR1" is bleaching and recovery 2014-2015
  filter(Period=="2014 Oct" | Period=="2015 Feb") 

# lilipuna bleaching-recovery2
Manova.Lil.BR2.lev.dat<-Manova.Lil.lev.dat %>% # "BR2" is bleaching and recovery 2014-2015
  filter(Period=="2015 Oct" | Period=="2016 Feb") 

Manova.Lil.BR1.lev.dat$Period<-droplevels(Manova.Lil.BR1.lev.dat$Period)
Manova.Lil.BR1.lev.dat$Site<-droplevels(Manova.Lil.BR1.lev.dat$Site)
Manova.Lil.BR1.lev.dat$Period.dom<-droplevels(Manova.Lil.BR1.lev.dat$Period.dom)

Manova.Lil.BR2.lev.dat$Period<-droplevels(Manova.Lil.BR2.lev.dat$Period)
Manova.Lil.BR2.lev.dat$Site<-droplevels(Manova.Lil.BR2.lev.dat$Site)
Manova.Lil.BR2.lev.dat$Period.dom<-droplevels(Manova.Lil.BR2.lev.dat$Period.dom)

# BR1 just data (sans factors)
Manova.Lil.BR1.dat<-Manova.Lil.BR1.lev.dat %>% 
  select(-c(Period, Status, Site, dom, Period.dom)) #  sans factors

# BR2 just data (sans factors)
Manova.Lil.BR2.dat<-Manova.Lil.BR2.lev.dat %>% 
  select(-c(Period, Status, Site, dom, Period.dom)) #  sans factors


############################
############### 
##  Run BR1 2014-2015 MANOVA Lilipuna

## PERMANOVA
Manova.Lil.BR1.dat.scaled <- scale(Manova.Lil.BR1.dat, center = TRUE, scale = TRUE) 

MVA.Lilipuna.BR1<-adonis2(Manova.Lil.BR1.dat.scaled~Period*dom, data=Manova.Lil.BR1.lev.dat,
                          permutations=1000, method="euclidian")

## NMDS plots
NMDS.Lil.BR1<-metaMDS(Manova.Lil.BR1.dat, distane='bray', k=2, trymax=100) 
stressplot(NMDS.Lil.BR1)
factor.df.Lil.BR1<-Manova.Lil.BR1.lev.dat %>% 
  select(c(Period, Site, Status, dom, Period.dom)) # get just the factors

NMDS.df.Lil.BR1<-data.frame(x=NMDS.Lil.BR1$point[,1], y=NMDS.Lil.BR1$point[,2], 
                    Period=as.factor(factor.df.Lil.BR1[,1]),
                    Site=as.factor(factor.df.Lil.BR1[,2]),
                    Status=as.factor(factor.df.Lil.BR1[,3]),
                    dom=as.factor(factor.df.Lil.BR1[,4]),
                    Period.dom=as.factor(factor.df.Lil.BR1[,5]))

colnames(NMDS.df.Lil.BR1)[1:2]<-c("MDS1", "MDS2")

# vectors for variables
vars.Lil.BR1<-envfit(NMDS.Lil.BR1, Manova.Lil.BR1.dat, permu=999)
#scores(vars, "vectors")



#########################
########  Run 2015-2016 MANOVA Lilipuna
## PERMANOVA
Manova.Lil.BR2.dat.scaled <- scale(Manova.Lil.BR2.dat, center = TRUE, scale = TRUE) 

MVA.Lilipuna.BR2<-adonis2(Manova.Lil.BR2.dat.scaled~Period*dom, data=Manova.Lil.BR2.lev.dat,
                          permutations=1000, method="euclidian")

## NMDS
NMDS.Lil.BR2<-metaMDS(Manova.Lil.BR2.dat, distane='bray', k=2, trymax=100) 
stressplot(NMDS.Lil.BR2)
factor.df.Lil.BR2<-Manova.Lil.BR2.lev.dat %>% 
  select(c(Period, Site, Status, dom, Period.dom)) # get just the factors

NMDS.df.Lil.BR2<-data.frame(x=NMDS.Lil.BR2$point[,1], y=NMDS.Lil.BR2$point[,2], 
                    Period=as.factor(factor.df.Lil.BR2[,1]),
                    Site=as.factor(factor.df.Lil.BR2[,2]),
                    Status=as.factor(factor.df.Lil.BR2[,3]),
                    dom=as.factor(factor.df.Lil.BR2[,4]),
                    Period.dom=as.factor(factor.df.Lil.BR2[,5]))

colnames(NMDS.df.Lil.BR2)[1:2]<-c("MDS1", "MDS2")

# vectors for variables
vars.Lil.BR2<-envfit(NMDS.Lil.BR2, Manova.Lil.BR2.dat, permu=999)
#scores(vars, "vectors")


########## ########## ########## ########## 
########## Reef 14 2014-2015
## Reef 14 only #######

## all data with Reef14 only
Manova.R14.lev.dat<-df.MV.noNA %>% # "lev.dat" has levels and data
  filter(Site=="Reef 14") 

# Reef 14 bleaching-recovery1
Manova.R14.BR1.lev.dat<-Manova.R14.lev.dat %>% # "BR1" is bleachigng and recovery 2014-2015
  filter(Period=="2014 Oct" | Period=="2015 Feb") 

# Reef 14 bleaching-recovery2
Manova.R14.BR2.lev.dat<-Manova.R14.lev.dat %>% # "BR2" is bleachigng and recovery 2014-2015
  filter(Period=="2015 Oct" | Period=="2016 Feb") 

Manova.R14.BR1.lev.dat$Period<-droplevels(Manova.R14.BR1.lev.dat$Period)
Manova.R14.BR1.lev.dat$Site<-droplevels(Manova.R14.BR1.lev.dat$Site)
Manova.R14.BR1.lev.dat$Period.dom<-droplevels(Manova.R14.BR1.lev.dat$Period.dom)

Manova.R14.BR2.lev.dat$Period<-droplevels(Manova.R14.BR2.lev.dat$Period)
Manova.R14.BR2.lev.dat$Site<-droplevels(Manova.R14.BR2.lev.dat$Site)
Manova.R14.BR2.lev.dat$Period.dom<-droplevels(Manova.R14.BR2.lev.dat$Period.dom)

# BR1 just data (sans factors)
Manova.R14.BR1.dat<-Manova.R14.BR1.lev.dat %>% 
  select(-c(Period, Status, Site, dom, Period.dom)) #  sans factors

# BR2 just data (sans factors)
Manova.R14.BR2.dat<-Manova.R14.BR2.lev.dat %>% 
  select(-c(Period, Status, Site, dom, Period.dom)) #  sans factors


############################# 
############### 
##  Run BR1 2014-2015 MANOVA Reef 14
Manova.R14.BR1.dat.scaled <- scale(Manova.R14.BR1.dat, center = TRUE, scale = TRUE) 

## PERMANOVA
MVA.R14.BR1<-adonis2(Manova.R14.BR1.dat.scaled~Period*dom, data=Manova.R14.BR1.lev.dat,permutations=1000, 
                method="euclidian")

## NMDS
NMDS.R14.BR1<-metaMDS(Manova.R14.BR1.dat, distane='bray', k=2, trymax=100) 
stressplot(NMDS.R14.BR1)
factor.df.R14.BR1<-Manova.R14.BR1.lev.dat %>% 
  select(c(Period, Site, Status, dom, Period.dom)) # get just the factors

NMDS.df.R14.BR1<-data.frame(x=NMDS.R14.BR1$point[,1], y=NMDS.R14.BR1$point[,2], 
                    Period=as.factor(factor.df.R14.BR1[,1]),
                    Site=as.factor(factor.df.R14.BR1[,2]),
                    Status=as.factor(factor.df.R14.BR1[,3]),
                    dom=as.factor(factor.df.R14.BR1[,4]),
                    Period.dom=as.factor(factor.df.R14.BR1[,5]))

colnames(NMDS.df.R14.BR1)[1:2]<-c("MDS1", "MDS2")

# vectors for variables
vars.R14.BR1<-envfit(NMDS.R14.BR1, Manova.R14.BR1.dat, permu=999)
#scores(vars, "vectors")



############################# 
########  Run 2015-2016 MANOVA R14
Manova.R14.BR2.dat.scaled <- scale(Manova.R14.BR2.dat, center = TRUE, scale = TRUE) 

## PERMANOVA
MVA.R14.BR2<-adonis2(Manova.R14.BR2.dat.scaled~Period*dom, data=Manova.R14.BR2.lev.dat,permutations=1000, 
                method="euclidean")

## NMDS
NMDS.R14.BR2<-metaMDS(Manova.R14.BR2.dat, distane='bray', k=2, trymax=100) 
stressplot(NMDS.R14.BR2)
factor.df.R14.BR2<-Manova.R14.BR2.lev.dat %>% 
  select(c(Period, Site, Status, dom, Period.dom)) # get just the factors

NMDS.df.R14.BR2<-data.frame(x=NMDS.R14.BR2$point[,1], y=NMDS.R14.BR2$point[,2], 
                    Period=as.factor(factor.df.R14.BR2[,1]),
                    Site=as.factor(factor.df.R14.BR2[,2]),
                    Status=as.factor(factor.df.R14.BR2[,3]),
                    dom=as.factor(factor.df.R14.BR2[,4]),
                    Period.dom=as.factor(factor.df.R14.BR2[,5]))

colnames(NMDS.df.R14.BR2)[1:2]<-c("MDS1", "MDS2")

# vectors for variables
vars.R14.BR2<-envfit(NMDS.R14.BR2, Manova.R14.BR2.dat, permu=999)
#scores(vars, "vectors")

code for NMDS figures separated by sites in each time period.

# for exporting 4 panel plot
par(mfcol=c(2,2), mar=c(4,4,1,1))

###################################### plots
###################################### start with Lilipuna 2014-2015, then 2015-2016... then Reef 14

############### Lilipuna 2014-2015

NMDS.df.Lil.BR1$Period.dom<-factor(NMDS.df.Lil.BR1$Period.dom, levels =
                                     c("2014 Oct.C", "2014 Oct.D", "2015 Feb.C", "2015 Feb.D"))
groups<-NMDS.df.Lil.BR1$Period.dom

MDS.lev.leg<-c("C-Bleach '14", "D-Bleach '14", "C-Recov '15", "D-Recov '15")
MDS.colors<-c("lightpink", "slategray1", "firebrick3", "dodgerblue3")

### 
NMDS.plot<-ordiplot(NMDS.Lil.BR1, type="n", main=substitute(paste("Lilipuna")), cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.4, 0.4), xlim=c(-0.6, 0.5), xaxt="n", xlab="")
axis(side = 1, labels = FALSE, tck = -0.05)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", cex=0.8, pch=16, col=MDS.colors[groups])
ordihull(NMDS.plot, groups=NMDS.df.Lil.BR1$Period.dom, draw="polygon", alpha=100, col=MDS.colors, border=MDS.colors)
legend("topright", legend=MDS.lev.leg, inset=c(0.7, 0), x.intersp=0.9, y.intersp=1, cex=0.7, pch=16, col=MDS.colors, pt.cex=1, bty="n")
par.new=T
plot(vars.Lil.BR1, col="black", p.max=0.05, cex=0.8, lwd=1)


############### Lilipuna 2015-2016

NMDS.df.Lil.BR2$Period.dom<-factor(NMDS.df.Lil.BR2$Period.dom, levels =
                                     c("2015 Oct.C", "2015 Oct.D", "2016 Feb.C", "2016 Feb.D"))
groups<-NMDS.df.Lil.BR2$Period.dom

MDS.lev.leg<-c("C-Bleach '15", "D-Bleach '15", "C-Recov '16", "D-Recov '16")
MDS.colors<-c("lightpink", "slategray1", "firebrick3", "dodgerblue3")

NMDS.plot<-ordiplot(NMDS.Lil.BR2, type="n", cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.4, 0.4), xlim=c(-0.6, 0.5))
axis(side = 2, labels = FALSE, tck = -0.05)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", cex=0.8, pch=16,  col=MDS.colors[groups])
ordihull(NMDS.plot, groups=NMDS.df.Lil.BR2$Period.dom, draw="polygon", alpha=100, col=MDS.colors, border=MDS.colors)
legend("topright", legend=MDS.lev.leg, inset=c(0.7, 0), x.intersp=0.9, y.intersp=1, cex=0.7, pch=16, col=MDS.colors, pt.cex=1, bty="n")
par.new=T
plot(vars.Lil.BR2, col="black", p.max=0.05, cex=0.8, lwd=1)


############### ############### ############### ############### ############### 
############### Reef 14 2014-2015

NMDS.df.R14.BR1$Period.dom<-factor(NMDS.df.R14.BR1$Period.dom, levels =
                                     c("2014 Oct.C", "2014 Oct.D", "2015 Feb.C", "2015 Feb.D"))
groups<-NMDS.df.R14.BR1$Period.dom

MDS.lev.leg<-c("C-Bleach '14", "D-Bleach '14", "C-Recov '15", "D-Recov '15")
MDS.colors<-c("lightpink", "slategray1", "firebrick2", "dodgerblue2")

NMDS.plot<-ordiplot(NMDS.R14.BR1, type="n", main=substitute(paste("Reef 14")), cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.4, 0.4), xlim=c(-0.6, 0.5), yaxt="n", ylab="", xaxt="n", xlab="")
axis(side = 1, labels = FALSE, tck = -0.05)
axis(side = 2, labels = FALSE, tck = -0.05)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", cex=0.8, pch=16,  col=MDS.colors[groups])
ordihull(NMDS.plot, groups=NMDS.df.R14.BR1$Period.dom, draw="polygon", alpha=100, col=MDS.colors, border=MDS.colors)
par.new=T
plot(vars.R14.BR1, col="black", p.max=0.05, cex=0.8, lwd=1)


############### Reef 14 2015-2016

NMDS.df.R14.BR2$Period.dom<-factor(NMDS.df.R14.BR2$Period.dom, levels =
                                     c("2015 Oct.C", "2015 Oct.D", "2016 Feb.C", "2016 Feb.D"))
groups<-NMDS.df.R14.BR2$Period.dom

MDS.lev.leg<-c("C-Bleach '15", "D-Bleach '15", "C-Recov '16", "D-Recov '16")
MDS.colors<-c("lightpink", "slategray1", "firebrick2", "dodgerblue2")

NMDS.plot<-ordiplot(NMDS.R14.BR2, type="n", cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.4, 0.4), xlim=c(-0.6, 0.5), yaxt="n", ylab="")
axis(side = 2, labels = FALSE, tck = -0.05)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", cex=0.8, pch=16,  col=MDS.colors[groups])
ordihull(NMDS.plot, groups=NMDS.df.R14.BR2$Period.dom, draw="polygon", alpha=100, col=MDS.colors, border=MDS.colors)
par.new=T
plot(vars.R14.BR2, col="black", p.max=0.05, cex=0.8, lwd=1)
**Figure 4.** Non-metric multidimensional scaling (NMDS) analyses identify coral physiotypes separated by dominant symbiont type (colors), bleaching-recovery periods (lighter and darker shades), sites (columns), and events (rows).

Figure 4. Non-metric multidimensional scaling (NMDS) analyses identify coral physiotypes separated by dominant symbiont type (colors), bleaching-recovery periods (lighter and darker shades), sites (columns), and events (rows).

dev.copy(pdf, "figures/4panel_NMDS.pdf", height=6, width=7, useDingbats=FALSE)
dev.off() 

Univariate analyses

Figures

# side by side figures separated by Sites 
### Physiological figures

color.scheme<-c("gray85", "gray65", "gray50", "gray20", 
                "firebrick1", #C Lilipuna
                "firebrick3", #C Reef 14
                "steelblue1", # D Lilipuna
                "steelblue3") # D Reef 14

break.points<-c("C.Lilipuna.Pre", "C.Reef 14.Pre", "D.Lilipuna.Pre", 
                "D.Reef 14.Pre", "C.Lilipuna", "C.Reef 14", "D.Lilipuna", "D.Reef 14")
legend.names<-c("C-Lilipuna Pre", "C-Reef 14 Pre", "D-Lilipuna Pre", "D-Reef 14 Pre",
                "C-Lilipuna", "C-Reef 14", "D-Lilipuna", "D-Reef 14")
axis.labels<-c("Prev", "Oct'14 \nBl","Feb'15 \nRc","Oct'15 \nBl", "Feb'16 \nRc")


Fig.formatting<-(theme_classic()) +
  theme(text=element_text(size=6),
  axis.line=element_blank(),
  legend.text.align = 0,
  legend.text=element_text(size=5),
    #legend.title = element_blank(),
      panel.border = element_rect(fill=NA, colour = "black", size=0.7),
      aspect.ratio=1, 
    axis.ticks.length=unit(0.2, "cm"),
    axis.text.y=element_text(
      margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=6), 
    axis.text.x=element_text(
      margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=6)) +
  theme(legend.key.size = unit(0.4, "cm")) +
  theme(aspect.ratio=1) +
  theme(panel.spacing=unit(c(0, 0, 0, 0), "cm"))


pd <- position_dodge(0.15) #offset for error bars


########################
#### graph as category of C/D
######################## 

fig.df.phys$int<-factor(fig.df.phys$int, levels=c("C.Lilipuna.Pre", "C.Reef 14.Pre", "D.Lilipuna.Pre", "D.Reef 14.Pre", "C.Lilipuna", "C.Reef 14", "D.Lilipuna", "D.Reef 14"))


########################
########################
########################
# Code for generating individual figures


###################
# zoox 
zoox.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=zoox/10^6, group=int, color=int)) + geom_errorbar(aes(ymin=zoox/10^6-zoox.SE/10^6, ymax=zoox/10^6+zoox.SE/10^6), size=.4, width=0, position=pd) + 
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5, pch=19) +
  coord_cartesian(ylim=c(0, 3)) +
  ylab(expression(paste(~("symbionts") ~(x10^6~cells~cm^-2), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Symbiont:Site",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting



########## chla cm2 ############
chla.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=chla, group=int, color=int)) + geom_errorbar(aes(ymin=chla-chla.SE, ymax=chla+chla.SE), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(0, 8)) +
  ylab(expression(paste("Chlorophyll", ~italic("a"), ~(mu*g~cm^-2), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Symbiont:Site",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting


########## chla cell ############
chlacell.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=chlacell, group=int, color=int)) + geom_errorbar(aes(ymin=chlacell-chlacell.SE, ymax=chlacell+chlacell.SE), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(0, 8)) +
  ylab(expression(paste("Chlorophyll", ~italic("a"), ~(pg~cell^-1), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Symbiont:Site",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting



########## protein ############
prot.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=prot, group=int, color=int)) + geom_errorbar(aes(ymin=prot-prot.SE, ymax=prot+prot.SE), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(0.2, 0.7)) +
  ylab(expression(paste("Protein", ~(mg~cm^-2), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Symbiont:Site",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting


########## AFDW ############
AFDW.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=AFDW, group=int, color=int)) + geom_errorbar(aes(ymin=AFDW-AFDW.SE, ymax=AFDW+AFDW.SE), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(10, 60)) +
  ylab(expression(paste("Total biomass", ~(mg~cm^-2), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Symbiont:Site",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting

########################################################################
########################################################################


### ### ### ### ### Immunity figures ### ### ### ### ### 

color.scheme.immuno<-c("gray70", "gray50", 
                "firebrick1", #C Lilipuna
                "firebrick3", #C Reef 14
                "steelblue1", # D Lilipuna
                "steelblue3") # D Reef 14


break.points.immuno<-c("unident.Lilipuna.Pre", "unident.Reef 14.Pre", 
                       "C.Lilipuna", "C.Reef 14", "D.Lilipuna", "D.Reef 14")
legend.names.immuno<-c("unk-Lilipuna Pre", "unk-Reef 14 Pre",
                       "C-Lilipuna", "C-Reef 14", "D-Lilipuna", "D-Reef 14")
axis.labels.immuno<-c("Prev", "Oct'14 \nBl","Feb'15 \nRc","Oct'15 \nBl", "Feb'16 \nRc")


########## Prophenoloxidase (PPO) ############
PPO.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=PPO, group=int, color=int)) + geom_errorbar(aes(ymin=PPO-PPO.SE, ymax=PPO+PPO.SE), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(0, 2)) +
  ylab(expression(paste("PPO Activity", ~(Delta~Abs~"490nm"~min^-1~mg~prot^-1), sep=""))) +
  scale_color_manual(values=c(color.scheme.immuno),
                     name= "Symbiont:Site",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting


### Melanin ###
MEL.lab<-aggregate(MEL~Site+Status+Period,data=data, mean)
MEL.labSE<-aggregate(MEL~Site+Status+Period,data=data, std.error)
MEL.df<-cbind(MEL.lab, MEL.labSE[c(4,0)]); colnames(MEL.df[,4:5])<-c("MEL", "MEL.labSE")
MEL.df$Period<-factor(MEL.df$Period, levels=c("2014 Field.Feb", "2014 Lab", "2014 Oct", "2015 Feb", "2015 Oct", "2016 Feb"))


######## Melanin (MEL) #############
MEL.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=MEL, group=int, color=int)) + geom_errorbar(aes(ymin=MEL-MEL.SE, ymax=MEL+MEL.SE), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(0, 0.02)) +
  ylab(expression(paste("MEL", ~(Abs~"490nm"~mg~tissue^-1), sep="")))+
  scale_color_manual(values=c(color.scheme.immuno),
                     name= "Symbiont:Site",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting


########## Peroxidase (POX) ##############
POX.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=POX, group=int, color=int)) + geom_errorbar(aes(ymin=POX-POX.SE, ymax=POX+POX.SE), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(0, 1.25)) +
  ylab(expression(paste("POX activity", ~(Delta~Abs~"470nm"~min^-1~mg~prot^-1), sep=""))) +
  scale_color_manual(values=c(color.scheme.immuno),
                     name= "Symbiont:Site",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting


###### Catalase (CAT) ########
CAT.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=CAT, group=int, color=int)) + geom_errorbar(aes(ymin=CAT-CAT.SE, ymax=CAT+CAT.SE), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(0, 500)) +
  ylab(expression(paste("CAT activity", ~(min^-1~mg~prot^-1), sep="")))+
  scale_color_manual(values=c(color.scheme.immuno),
                     name= "Symbiont:Site",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting


######### Superoxide dismutase (SOD) ########
SOD.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=SOD/10^3, group=int, color=int)) + geom_errorbar(aes(ymin=SOD/10^3-SOD.SE/10^3, ymax=SOD/10^3+SOD.SE/10^3), size=.4, width=0, position=pd) +
  geom_line(aes(linetype=SimpleSite), position=pd, lwd=0.4) +
  scale_linetype_manual(values=c("solid", "longdash"))+
  geom_point(position=pd, size=1.5) +
  coord_cartesian(ylim=c(10, 40)) +
  ylab(expression(paste("SOD activity", ~(x10^3~mg~prot^-1), sep=""))) +
  scale_color_manual(values=c(color.scheme.immuno),
                     name= "Symbiont:Site",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Period", 
      labels=axis.labels)+
  Fig.formatting

The combined plots for physiology and immunity data.

# legend plot
fig.legend.phys<-ggplot(data=fig.df.phys, aes(x=Period, y=zoox, group=int, color=int, linetype=int)) +
  geom_line() +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 1)) +
  scale_color_manual(values=c(color.scheme),
                     name= "Symbiont:Site",
                     breaks=break.points,
                     labels=legend.names) +
  scale_linetype_manual(values=c(0,0,0,0,1,3,1,3),
                     name= "Symbiont:Site",
                     breaks=break.points,
                     labels=legend.names) +
  theme_void()

Phys.figs<-plot_grid( (zoox.fig + theme(legend.position = "none")), 
           (chla.fig + theme(legend.position = "none")),
           (chlacell.fig + theme(legend.position = "none")), 
           (prot.fig+ theme(legend.position = "none")),
           (AFDW.fig+ theme(legend.position = "none")),
           fig.legend.phys,
          labels = c("a", "b", "c", "d", "e",""), ncol = 3, label_size=8)
Phys.figs
**Supplementary Figure 2.** Physiological metrics for *M. capitata* corals dominated dominated by *Cladocopium* sp. or *Durusdinium* sp. symbionts (C or D) from two reefs in Kāne‘ohe Bay during repeated bleaching and recovery periods.

Supplementary Figure 2. Physiological metrics for M. capitata corals dominated dominated by Cladocopium sp. or Durusdinium sp. symbionts (C or D) from two reefs in Kāne‘ohe Bay during repeated bleaching and recovery periods.

dev.copy(pdf, "figures/Phys.figs.pdf", encod="MacRoman", height=5, width=8)
dev.off()
fig.legend.immun<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=SOD, group=int, color=int, linetype=int)) +
  geom_point(position=pd, size=3) +
  geom_line() +
  coord_cartesian(ylim=c(0, 1)) +
  scale_color_manual(values=c(color.scheme.immuno),
                     name= "Symbiont:Site",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_linetype_manual(values=c(0,0,1,3,1,3),
                     name= "Symbiont:Site",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  theme_void()


Immuno.figs<-plot_grid( (MEL.fig + theme(legend.position = "none")),
           (PPO.fig + theme(legend.position = "none")), 
           (POX.fig+ theme(legend.position = "none")),
           (CAT.fig+ theme(legend.position = "none")),
           SOD.fig + theme(legend.position = "none"),
           fig.legend.immun,
          labels = c("a", "b", "c", "d", "e",""), ncol = 3, label_size=8)
Immuno.figs
**Supplementary Figure 3.** Immunity metrics for *M. capitata* corals dominated dominated by *Cladocopium* sp. or *Durusdinium* sp. symbionts (C or D) from two reefs in Kāne‘ohe Bay during repeated bleaching and recovery periods.

Supplementary Figure 3. Immunity metrics for M. capitata corals dominated dominated by Cladocopium sp. or Durusdinium sp. symbionts (C or D) from two reefs in Kāne‘ohe Bay during repeated bleaching and recovery periods.

dev.copy(pdf, "figures/Immuno.figs.pdf", encod="MacRoman", height=5, width=8)
dev.off()

Models

Running models of response variables. First check for assumptions of ANOVA and apply transformations where necessary.

lm(Y~Period x Site x dom, data=model.data, na.action=na.exclude)
Period is the 4 Periods, Site is the two reefs, dom is the symbiont community (Cladooopium- or Durusdinium-dominated, ie., C- or D-dominated).

Physiology models

Run linear models for response variables

Symbionts cm-2
models and post-hoc comparisons

# symbionts ---- 
Y<-model.data$zoox
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                    Sum Sq  Df F value   Pr(>F)    
## Period          1.108e+13   3   9.470 5.37e-06 ***
## Site            3.068e+10   1   0.079  0.77928    
## dom             5.640e+13   1 144.647  < 2e-16 ***
## Period:Site     4.604e+12   3   3.936  0.00888 ** 
## Period:dom      3.461e+12   3   2.959  0.03260 *  
## Site:dom        8.347e+11   1   2.141  0.14447    
## Period:Site:dom 4.342e+12   3   3.712  0.01198 *  
## Residuals       1.174e+14 301                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
multcomp::cld(posthoc, Letters=letters)
##  Period    lsmean    SE  df lower.CL upper.CL .group
##  2014 Oct 1215332 71271 301  1075079  1355585  a    
##  2015 Oct 1351758 69900 301  1214204  1489313  ab   
##  2015 Feb 1491352 74355 301  1345031  1637674   bc  
##  2016 Feb 1745350 72651 301  1602383  1888318    c  
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom)
multcomp::cld(posthoc, Letters=letters)
##  dom  lsmean    SE  df lower.CL upper.CL .group
##  C   1024246 53718 301   918535  1129956  a    
##  D   1877651 48036 301  1783122  1972180   b   
## 
## Results are averaged over the levels of: Period, Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site      lsmean     SE  df lower.CL upper.CL .group
##  Reef 14  1074234 100820 301   875834  1272634  a    
##  Lilipuna 1356430 100766 301  1158136  1554724   b   
## 
## Period = 2015 Feb:
##  Site      lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna 1389494 111098 301  1170868  1608121  a    
##  Reef 14  1593210  98853 301  1398679  1787741  a    
## 
## Period = 2015 Oct:
##  Site      lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna 1216313  98853 301  1021782  1410844  a    
##  Reef 14  1487203  98853 301  1292672  1681735  a    
## 
## Period = 2016 Feb:
##  Site      lsmean     SE  df lower.CL upper.CL .group
##  Reef 14  1689582 100020 301  1492754  1886410  a    
##  Lilipuna 1801118 105396 301  1593712  2008525  a    
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom  lsmean     SE  df lower.CL upper.CL .group
##  C    731023 102581 301   529156   932891  a    
##  D   1699641  98971 301  1504877  1894404   b   
## 
## Period = 2015 Feb:
##  dom  lsmean     SE  df lower.CL upper.CL .group
##  C   1247815 118287 301  1015041  1480588  a    
##  D   1734890  90128 301  1557530  1912250   b   
## 
## Period = 2015 Oct:
##  dom  lsmean     SE  df lower.CL upper.CL .group
##  C    803195  96350 301   613590   992801  a    
##  D   1900322 101295 301  1700986  2099657   b   
## 
## Period = 2016 Feb:
##  dom  lsmean     SE  df lower.CL upper.CL .group
##  C   1314949 111229 301  1096064  1533834  a    
##  D   2175752  93491 301  1991774  2359730   b   
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site*dom|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site     dom  lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna C    634517 156105 301   327321   941714  a    
##  Reef 14  C    827529 133127 301   565551  1089507  ab   
##  Reef 14  D   1320939 151445 301  1022915  1618963   b   
##  Lilipuna D   2078343 127460 301  1827518  2329167    c  
## 
## Period = 2015 Feb:
##  Site     dom  lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna C   1214202 188270 301   843709  1584694  a    
##  Reef 14  C   1281428 143252 301   999525  1563330  a    
##  Lilipuna D   1564787 118005 301  1332569  1797006  ab   
##  Reef 14  D   1904993 136260 301  1636850  2173136   b   
## 
## Period = 2015 Oct:
##  Site     dom  lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna C    705542 136260 301   437399   973685  a    
##  Reef 14  C    900848 136260 301   632705  1168991  a    
##  Lilipuna D   1727084 143252 301  1445182  2008987   b   
##  Reef 14  D   2073559 143252 301  1791656  2355461   b   
## 
## Period = 2016 Feb:
##  Site     dom  lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna C   1297075 173183 301   956271  1637878  a    
##  Reef 14  C   1332823 139625 301  1058058  1607587  a    
##  Reef 14  D   2046342 143252 301  1764439  2328244   b   
##  Lilipuna D   2305162 120170 301  2068682  2541641   b   
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="zoox", par.strip.text=list(cex=0.7))

Chlorophyll a cm-2
models and post-hoc comparisons

# chla ---- 
Y<-model.data$chla
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period            85.3   3  11.986 1.96e-07 ***
## Site              23.9   1  10.049  0.00168 ** 
## dom                3.4   1   1.438  0.23137    
## Period:Site        7.0   3   0.988  0.39860    
## Period:dom        66.1   3   9.278 6.92e-06 ***
## Site:dom           0.1   1   0.029  0.86516    
## Period:Site:dom   12.7   3   1.790  0.14906    
## Residuals        714.4 301                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
multcomp::cld(posthoc, Letters=letters)
##  Period   lsmean    SE  df lower.CL upper.CL .group
##  2015 Oct   3.16 0.172 301     2.82     3.50  a    
##  2014 Oct   3.40 0.176 301     3.06     3.75  a    
##  2015 Feb   4.22 0.183 301     3.85     4.58   b   
##  2016 Feb   4.52 0.179 301     4.16     4.87   b   
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
multcomp::cld(posthoc, Letters=letters)
##  Site     lsmean    SE  df lower.CL upper.CL .group
##  Lilipuna   3.55 0.128 301     3.30     3.81  a    
##  Reef 14    4.09 0.123 301     3.85     4.34   b   
## 
## Results are averaged over the levels of: Period, dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom lsmean    SE  df lower.CL upper.CL .group
##  C     3.10 0.253 301     2.60     3.60  a    
##  D     3.71 0.244 301     3.23     4.19  a    
## 
## Period = 2015 Feb:
##  dom lsmean    SE  df lower.CL upper.CL .group
##  D     3.93 0.222 301     3.50     4.37  a    
##  C     4.50 0.292 301     3.92     5.07  a    
## 
## Period = 2015 Oct:
##  dom lsmean    SE  df lower.CL upper.CL .group
##  C     2.42 0.238 301     1.95     2.88  a    
##  D     3.90 0.250 301     3.41     4.39   b   
## 
## Period = 2016 Feb:
##  dom lsmean    SE  df lower.CL upper.CL .group
##  D     4.09 0.231 301     3.63     4.54  a    
##  C     4.94 0.274 301     4.40     5.48   b   
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="chla", par.strip.text=list(cex=0.7))

chlorophyll a cell-1
models and post-hoc comparisons

# chlacell ---- 
Y<-model.data$chlacell
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
##summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period           19.27   3   6.870 0.000173 ***
## Site              5.97   1   6.382 0.012044 *  
## dom             181.23   1 193.843  < 2e-16 ***
## Period:Site       4.06   3   1.449 0.228640    
## Period:dom       14.11   3   5.029 0.002051 ** 
## Site:dom          3.23   1   3.459 0.063903 .  
## Period:Site:dom   6.15   3   2.192 0.089017 .  
## Residuals       278.61 298                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
multcomp::cld(posthoc, Letters=letters)
##  Period   lsmean    SE  df lower.CL upper.CL .group
##  2015 Oct   2.80 0.110 298     2.58     3.01  a    
##  2016 Feb   2.87 0.114 298     2.65     3.10  a    
##  2015 Feb   2.96 0.115 298     2.74     3.19  a    
##  2014 Oct   3.48 0.110 298     3.27     3.70   b   
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
multcomp::cld(posthoc, Letters=letters)
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna   2.89 0.0814 298     2.73     3.05  a    
##  Reef 14    3.17 0.0774 298     3.01     3.32   b   
## 
## Results are averaged over the levels of: Period, dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom)
multcomp::cld(posthoc, Letters=letters)
##  dom lsmean     SE  df lower.CL upper.CL .group
##  D     2.25 0.0747 298     2.10     2.40  a    
##  C     3.81 0.0840 298     3.64     3.98   b   
## 
## Results are averaged over the levels of: Period, Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom lsmean    SE  df lower.CL upper.CL .group
##  D     2.41 0.153 298     2.11     2.71  a    
##  C     4.55 0.159 298     4.24     4.87   b   
## 
## Period = 2015 Feb:
##  dom lsmean    SE  df lower.CL upper.CL .group
##  D     2.32 0.140 298     2.04     2.59  a    
##  C     3.61 0.183 298     3.25     3.97   b   
## 
## Period = 2015 Oct:
##  dom lsmean    SE  df lower.CL upper.CL .group
##  D     2.28 0.159 298     1.96     2.59  a    
##  C     3.32 0.151 298     3.02     3.62   b   
## 
## Period = 2016 Feb:
##  dom lsmean    SE  df lower.CL upper.CL .group
##  D     1.99 0.145 298     1.70     2.27  a    
##  C     3.76 0.177 298     3.41     4.10   b   
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="AFDW", par.strip.text=list(cex=0.7))

Protein cm-2
models and post-hoc comparisons

# prot ---- 
Y<-model.data$prot
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value  Pr(>F)   
## Period           0.058   3   0.722 0.53971   
## Site             0.053   1   1.980 0.16042   
## dom              0.126   1   4.706 0.03085 * 
## Period:Site      0.412   3   5.135 0.00178 **
## Period:dom       0.035   3   0.440 0.72477   
## Site:dom         0.011   1   0.394 0.53082   
## Period:Site:dom  0.048   3   0.596 0.61820   
## Residuals        8.020 300                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~dom)
multcomp::cld(posthoc, Letters=letters)
##  dom lsmean     SE  df lower.CL upper.CL .group
##  C    0.441 0.0141 300    0.413    0.469  a    
##  D    0.480 0.0126 300    0.456    0.505   b   
## 
## Results are averaged over the levels of: Period, Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Reef 14   0.408 0.0269 300    0.355    0.461  a    
##  Lilipuna  0.539 0.0264 300    0.487    0.591   b   
## 
## Period = 2015 Feb:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna  0.414 0.0291 300    0.356    0.471  a    
##  Reef 14   0.487 0.0259 300    0.436    0.538  a    
## 
## Period = 2015 Oct:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Reef 14   0.437 0.0259 300    0.386    0.488  a    
##  Lilipuna  0.459 0.0259 300    0.408    0.510  a    
## 
## Period = 2016 Feb:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Reef 14   0.455 0.0262 300    0.404    0.507  a    
##  Lilipuna  0.486 0.0276 300    0.432    0.540  a    
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="prot", par.strip.text=list(cex=0.7))

Total biomass cm-2
models and post-hoc comparisons

# prot ---- 
Y<-model.data$AFDW
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period           7.363   3  78.602  < 2e-16 ***
## Site             0.603   1  19.312 1.54e-05 ***
## dom              0.355   1  11.359 0.000849 ***
## Period:Site      1.213   3  12.954 5.57e-08 ***
## Period:dom       0.346   3   3.696 0.012234 *  
## Site:dom         0.085   1   2.738 0.099031 .  
## Period:Site:dom  0.065   3   0.698 0.554063    
## Residuals        9.398 301                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Site|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Reef 14    2.05 0.0285 301     1.99     2.10  a    
##  Lilipuna   2.20 0.0285 301     2.14     2.26   b   
## 
## Period = 2015 Feb:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Lilipuna   1.89 0.0314 301     1.83     1.95  a    
##  Reef 14    2.01 0.0280 301     1.96     2.07   b   
## 
## Period = 2015 Oct:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Reef 14    2.31 0.0280 301     2.25     2.36  a    
##  Lilipuna   2.44 0.0280 301     2.39     2.50   b   
## 
## Period = 2016 Feb:
##  Site     lsmean     SE  df lower.CL upper.CL .group
##  Reef 14    2.11 0.0283 301     2.06     2.17  a    
##  Lilipuna   2.30 0.0298 301     2.25     2.36   b   
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom lsmean     SE  df lower.CL upper.CL .group
##  C     2.11 0.0290 301     2.05     2.17  a    
##  D     2.14 0.0280 301     2.08     2.19  a    
## 
## Period = 2015 Feb:
##  dom lsmean     SE  df lower.CL upper.CL .group
##  C     1.94 0.0335 301     1.87     2.01  a    
##  D     1.96 0.0255 301     1.91     2.01  a    
## 
## Period = 2015 Oct:
##  dom lsmean     SE  df lower.CL upper.CL .group
##  C     2.28 0.0273 301     2.23     2.34  a    
##  D     2.47 0.0287 301     2.41     2.52   b   
## 
## Period = 2016 Feb:
##  dom lsmean     SE  df lower.CL upper.CL .group
##  C     2.19 0.0315 301     2.13     2.26  a    
##  D     2.23 0.0265 301     2.17     2.28  a    
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="AFDW", par.strip.text=list(cex=0.7))

Immunity models

Melanin (MEL)
models and post-hoc comparisons

# MEL ---- 
Y<-model.data$MEL
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df  F value   Pr(>F)    
## Period          1.7128   3 1133.636  < 2e-16 ***
## Site            0.0006   1    1.112    0.292    
## dom             0.0000   1    0.000    0.987    
## Period:Site     0.0155   3   10.241 1.95e-06 ***
## Period:dom      0.0019   3    1.276    0.283    
## Site:dom        0.0001   1    0.198    0.657    
## Period:Site:dom 0.0003   3    0.170    0.917    
## Residuals       0.1496 297                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
multcomp::cld(posthoc, Letters=letters)
##  Period   lsmean      SE  df lower.CL upper.CL .group
##  2015 Feb  0.137 0.00267 297    0.132    0.142  a    
##  2016 Feb  0.216 0.00261 297    0.210    0.221   b   
##  2015 Oct  0.230 0.00256 297    0.225    0.235    c  
##  2014 Oct  0.345 0.00257 297    0.340    0.350     d 
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site     lsmean      SE  df lower.CL upper.CL .group
##  Reef 14   0.344 0.00362 297    0.337    0.351  a    
##  Lilipuna  0.345 0.00365 297    0.338    0.353  a    
## 
## Period = 2015 Feb:
##  Site     lsmean      SE  df lower.CL upper.CL .group
##  Lilipuna  0.126 0.00399 297    0.118    0.134  a    
##  Reef 14   0.149 0.00355 297    0.142    0.156   b   
## 
## Period = 2015 Oct:
##  Site     lsmean      SE  df lower.CL upper.CL .group
##  Reef 14   0.222 0.00360 297    0.215    0.229  a    
##  Lilipuna  0.239 0.00365 297    0.232    0.246   b   
## 
## Period = 2016 Feb:
##  Site     lsmean      SE  df lower.CL upper.CL .group
##  Lilipuna  0.212 0.00379 297    0.204    0.219  a    
##  Reef 14   0.219 0.00359 297    0.212    0.226  a    
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="prot", par.strip.text=list(cex=0.7))

Prophenoloxidase (PPO)
models and post-hoc comparisons

# PPO ---- 
Y<-model.data$PPO
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value Pr(>F)    
## Period           8.112   3 207.503 <2e-16 ***
## Site             0.055   1   4.227 0.0407 *  
## dom              0.054   1   4.135 0.0429 *  
## Period:Site      0.002   3   0.051 0.9848    
## Period:dom       0.020   3   0.510 0.6757    
## Site:dom         0.000   1   0.005 0.9419    
## Period:Site:dom  0.001   3   0.031 0.9927    
## Residuals        3.857 296                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
multcomp::cld(posthoc, Letters=letters)
##  Period   lsmean     SE  df lower.CL upper.CL .group
##  2015 Oct  0.652 0.0130 296    0.627    0.678  a    
##  2014 Oct  0.723 0.0132 296    0.698    0.749   b   
##  2016 Feb  0.759 0.0134 296    0.733    0.785   b   
##  2015 Feb  1.072 0.0136 296    1.045    1.099    c  
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
multcomp::cld(posthoc, Letters=letters)
##  Site     lsmean      SE  df lower.CL upper.CL .group
##  Lilipuna  0.788 0.00963 296    0.769    0.807  a    
##  Reef 14   0.815 0.00914 296    0.797    0.833   b   
## 
## Results are averaged over the levels of: Period, dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom)
multcomp::cld(posthoc, Letters=letters)
##  dom lsmean      SE  df lower.CL upper.CL .group
##  C    0.788 0.00986 296    0.769    0.807  a    
##  D    0.815 0.00888 296    0.798    0.833   b   
## 
## Results are averaged over the levels of: Period, Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="PPO", par.strip.text=list(cex=0.7))

Peroxidase (POX)
models and post-hoc comparisons

# POX ---- 
Y<-model.data$POX
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period           1.657   3  16.504 6.51e-10 ***
## Site             0.121   1   3.619   0.0581 .  
## dom              0.001   1   0.042   0.8382    
## Period:Site      0.089   3   0.884   0.4500    
## Period:dom       0.363   3   3.612   0.0138 *  
## Site:dom         0.018   1   0.547   0.4602    
## Period:Site:dom  0.014   3   0.138   0.9372    
## Residuals        9.502 284                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
multcomp::cld(posthoc, Letters=letters)
##  Period   lsmean     SE  df lower.CL upper.CL .group
##  2015 Oct  0.703 0.0224 284    0.659    0.747  a    
##  2014 Oct  0.835 0.0215 284    0.793    0.878   b   
##  2015 Feb  0.843 0.0219 284    0.800    0.886   b   
##  2016 Feb  0.912 0.0213 284    0.870    0.954   b   
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom lsmean     SE  df lower.CL upper.CL .group
##  D    0.804 0.0301 284    0.745    0.864  a    
##  C    0.866 0.0306 284    0.806    0.926  a    
## 
## Period = 2015 Feb:
##  dom lsmean     SE  df lower.CL upper.CL .group
##  D    0.830 0.0268 284    0.778    0.883  a    
##  C    0.855 0.0346 284    0.787    0.923  a    
## 
## Period = 2015 Oct:
##  dom lsmean     SE  df lower.CL upper.CL .group
##  C    0.642 0.0305 284    0.582    0.702  a    
##  D    0.765 0.0329 284    0.700    0.829   b   
## 
## Period = 2016 Feb:
##  dom lsmean     SE  df lower.CL upper.CL .group
##  D    0.890 0.0274 284    0.836    0.944  a    
##  C    0.933 0.0326 284    0.869    0.998  a    
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="POX", par.strip.text=list(cex=0.7))

Catalase (CAT)
models and post-hoc comparisons

######################## 
### immunology
########################

# CAT ---- 
Y<-model.data$CAT
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period          3048.5   3 117.049  < 2e-16 ***
## Site             185.2   1  21.334 5.80e-06 ***
## dom               81.0   1   9.328  0.00247 ** 
## Period:Site      228.3   3   8.766 1.39e-05 ***
## Period:dom        57.9   3   2.222  0.08568 .  
## Site:dom           1.6   1   0.183  0.66900    
## Period:Site:dom   22.5   3   0.862  0.46111    
## Residuals       2526.3 291                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
multcomp::cld(posthoc, Letters=letters)
##  Period   lsmean    SE  df lower.CL upper.CL .group
##  2015 Feb   9.35 0.351 291     8.66     10.0  a    
##  2016 Feb   9.36 0.343 291     8.69     10.0  a    
##  2014 Oct  12.42 0.339 291    11.75     13.1   b   
##  2015 Oct  17.18 0.349 291    16.50     17.9    c  
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
multcomp::cld(posthoc, Letters=letters)
##  Site     lsmean    SE  df lower.CL upper.CL .group
##  Reef 14    11.3 0.240 291     10.8     11.7  a    
##  Lilipuna   12.9 0.249 291     12.4     13.4   b   
## 
## Results are averaged over the levels of: Period, dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom)
multcomp::cld(posthoc, Letters=letters)
##  dom lsmean    SE  df lower.CL upper.CL .group
##  C     11.6 0.256 291     11.1     12.1  a    
##  D     12.6 0.232 291     12.1     13.0   b   
## 
## Results are averaged over the levels of: Period, Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
multcomp::cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site     lsmean    SE  df lower.CL upper.CL .group
##  Reef 14   11.93 0.481 291    10.99    12.88  a    
##  Lilipuna  12.91 0.480 291    11.97    13.86  a    
## 
## Period = 2015 Feb:
##  Site     lsmean    SE  df lower.CL upper.CL .group
##  Lilipuna   9.16 0.524 291     8.13    10.19  a    
##  Reef 14    9.54 0.466 291     8.62    10.45  a    
## 
## Period = 2015 Oct:
##  Site     lsmean    SE  df lower.CL upper.CL .group
##  Reef 14   14.97 0.500 291    13.98    15.95  a    
##  Lilipuna  19.40 0.486 291    18.44    20.35   b   
## 
## Period = 2016 Feb:
##  Site     lsmean    SE  df lower.CL upper.CL .group
##  Reef 14    8.56 0.472 291     7.64     9.49  a    
##  Lilipuna  10.16 0.497 291     9.18    11.14   b   
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="CAT", par.strip.text=list(cex=0.7))

Superoxide dismutase (SOD)
models and post-hoc comparisons

# SOD ---- 
Y<-model.data$SOD
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                    Sum Sq  Df F value Pr(>F)    
## Period          1.349e+10   3  83.207 <2e-16 ***
## Site            2.631e+08   1   4.867 0.0281 *  
## dom             8.110e+07   1   1.500 0.2216    
## Period:Site     9.381e+07   3   0.578 0.6296    
## Period:dom      2.319e+08   3   1.430 0.2340    
## Site:dom        2.041e+07   1   0.378 0.5394    
## Period:Site:dom 1.021e+08   3   0.630 0.5964    
## Residuals       1.616e+10 299                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
multcomp::cld(posthoc, Letters=letters)
##  Period   lsmean  SE  df lower.CL upper.CL .group
##  2014 Oct  15451 839 299    13799    17102  a    
##  2015 Feb  22587 875 299    20864    24310   b   
##  2015 Oct  29278 835 299    27635    30921    c  
##  2016 Feb  32585 855 299    30901    34268     d 
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
multcomp::cld(posthoc, Letters=letters)
##  Site     lsmean  SE  df lower.CL upper.CL .group
##  Reef 14   24030 589 299    22872    25189  a    
##  Lilipuna  25920 615 299    24710    27131   b   
## 
## Results are averaged over the levels of: Period, dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="SOD", par.strip.text=list(cex=0.7))

Schematic

Taking into account the observed patterns for bleaching and recovery in M. capitata we propose a shift in coral immunity and antioxidants in response to bleaching/recovery with a trend towards antioxidant reliance and increased constitutive responses.

Figure 5.. Schematic of coral responses to repeated bleaching and recovery.